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PREFACE 


The purpose of this study effort was to develop analytical models to describe the 
effects of residual accelerations on the experiments to be carried on the first U. 

S. Microgravity Lab mission (USML-1) and to test the accuracy of these models 
by comparing the pre-flight predicted effects with the post-flight measured effects. 
After surveying the experiments to be performed on USML-1, it became evident 
that the anticipated residual accelerations during the USML-1 mission were well 
below the threshold for most of the primary experiments and all of the secondary 
(Glovebox) experiments and that the only set of experiments that could provide 
quantifiable effects, and thus provide a definitive test of the analytical models, 
were the three melt growth experiments using the Bridgman-Stockbarger type 
Crystal Growth Furnace (CGF). This class of experiments is by far the most 
sensitive to low level quasi-steady accelerations that are unavoidable on space 
craft operating in low earth orbit. Because of this, they have been the drivers for 
the acceleration requirements imposed on the Space Station. Therefore, it is 
appropriate that the models on which these requirements are based are tested 
experimentally. Also, since solidification proceeds directionally over a long 
period of time, the solidified ingot provides a more or less continuous record of 
the effects from acceleration disturbances. 

The use of approximate analytical models is advantageous because of the insight 
they provide into the interaction of the many parameters involved in the 
redistribution of solute during the process. Numerical models are generally more 
accurate since they do not make as many simplifying assumptions, but they 
provide an answer for only a single point in parameter space and do not provide 
any indication of how the result might vary if one or more of the many input 
conditions were varied. Approximate analytical models provide scaling 
relationships that can cover a broad range of parameter space. If these scaling 
laws are verified by numerical computations at selected points to assure that the 
approximations used are still valid, they provide a very powerful tool for 
assessing the effects of residual accelerations on a variety of experiments 
without having to resort to a large number of numerical computations, and can 
provide a useful guide for experimenters wishing to minimize the unwanted 
effects from the residual accelerations. 

The first section of this report is devoted to the development of analytical models 
to describe the buoyancy-driven flows in closed containers with a variety of 
shapes and orientations relative to the acceleration vector. Primary emphasis 
was given to fully developed flows from quasi-steady accelerations, but methods 
for estimating the response times were also developed to extend these models to 
describe the effects of transient and periodic accelerations. A powerful integral 
theorem was derived which enabled the extension of simple one-dimensional 
models, for which closed-form analytical solutions exist, to two-dimensional 
models which have no known analytical solutions. The range of the two- 
dimensional models covers the case of low Rayleigh numbers, where the density 
field is decoupled from flow field, to higher Rayleigh numbers where such 
coupling must be considered. This allows the models to be extended nearly to 
unit gravity for the case of low Pr fluids, or to be applied to high Pr fluids as well 



as to solutal-d riven convection in a reduced gravity environment. The non-linear 
inertial terms were ignored, so the models are limited to flows where these terms 
are negligible. However, this is not a serious restriction for the purposes of this 
study. The effects of stabilizing and destabilizing density gradients were also 
considered along with the effects of applying magnetic fields to help suppress 
unwanted flows. These various models were tested extensively against 
numerical computations and were found to be accurate to within a factor of 2 and 
in most case accurate to within a few percent. 

Since USML-1 was to be a low-gravity emphasis mission, the orientation to the 
vehicle was chosen to orient the quasi-steady residual acceleration vector, 
primarily from atmospheric drag, along the furnace axis. Therefore, it was 
expected that the primary effect to be observed would be the radial segregation 
resulting from the radial thermal gradients in the furnace experiments. Therefore, 
the second section is devoted to modeling axisymmetrical flows and their effect 
on solute redistribution for a Bridgman growth system with low-level steady and 
time-dependent axial accelerations. This was accomplished by first using an 
approximate method to obtain the flow field and than using a perturbation method 
to get a first order correction to the concentration field. This method is applicable 
as long as the perturbation is small compared to the zero-order field. Again the 
results were tested against numerical calculations and were found to be 
reasonably accurate within the limits of the approximations used. 

As it turned out, there were some unanticipated accelerations during the mission 
that caused the acceleration vector at the furnace location to be more transverse 
than axial. This prompted an additional, more detailed study of the flows and 
their effects from transverse accelerations which is developed in Section 3. 

Since the Bridgman type experiments generally have a density profile that falls 
exponentially from the growth interface, emphasis was given to modeling the 
flows resulting from this type of density profile rather than using the simpler 
assumption of a linear profile which was the basis of most of the models 
developed in Section 1 . Also, to calculate the radial segregation at the growth 
interface, it is necessary to have an accurate description of the axial velocity 
profile. This was obtained using a novel approximation technique for the 
solution of the biharmonic equation that describes the stream function. A first 
order perturbation correction to the concentration field was then obtained as 
before. Again the results were verified by extensive comparison to numerical 
computations. 

Section 4 contains the pre-flight predictions of the expected effects of the 
anticipated accelerations. It should be remembered that since at the time these 
predictions were made it was expected that axial accelerations would be the 
dominant disturbance and that the study in Section 3 had not been undertaken. 
Therefore, these estimated are base primarily on the results of Section 2 and 
some crude estimates of the effects of transverse accelerations based on flow 
models in Section 1 . 

Section 5 contains more accurate after-the fact predictions of the effects of 
transverse accelerations using the results developed in Section 3. In comparing 
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th© two, it may b© sssn that th© crud© ©stimatss in Section 4 wore well within an 
order of magnitude of the more refined estimates. Unfortunately, experimental 
difficulties with the various experiments caused by the unanticipated transverse 
accelerations did not permit a definitive comparison of the models with 
experimental data. One of the lessons learned, however, was the extreme 
sensitivity of this type of experiment to quasi-steady transverse accelerations; 
especially for non-dilute alloy-type systems with high Schmidt numbers. Several 
possible methods such as density gradient stabilization and magnetic 
suppression were explored to help mitigate this sensitivity. 
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SECTION 1 


METHODS FOR ESTIMATING BUOYANCY-DRIVEN FLOWS 
IN CLOSED CAVITIES 


1.1 INTRODUCTION 


In many cases, it is desirable to be able to estimate the magnitude of natural or 
buoyancy-driven convective flows in a crystal growth or other solidification 
process in order to determine if heat and mass transport are dominated by 
convective or diffusive transport, to determine if flows will be steady or time- 
varying or to estimate the effectiveness of the strategy for controlling unwanted 
flows, whether by using strong magnetic fields or by conducting the experiment in 
reduced gravity. Such flows can, of course be computed by solving the full set of 
coupled Navier-Stokes equations in 2- or 3 dimensions, but such computations 
require a large computer, are expensive, and often give far more information 
than is really needed. Furthermore, such computations only give the solution to a 
specific problem and do not provide the insight necessary to judge the potential 
effects or varying the several parameters that determine the flow field 
Therefore, in order to develop a process, it is useful to have simple analytical 
models that are admittedly crude, but none-the-less can provide reasonably 
accurate predictions of the flows to be expected. 


Ostrach [1] pointed out the usefulness of applying scaling analysis to estimate 
flows in potential experiments to at least determine if the desired effects could be 
realized. Using this technique, he obtained 


v = Gr(v/L), Gr < 1 

v = Gr 1 ' 2 (v/L), Gr> 1 

v = (Gr/Pr) 1/2 (v/L), Gr>1 


and Pr<l; 

(i.i) 

and Pr<l; 

(1.2) 

and Pr > 1; 

(1.3) 


where v is the maximum velocity, Gr is the Grashof number (g Ap/p L 3 /v 2 ), Pr is 

the Prandtl number (v/k), v is the kinematic viscosity, k is the thermal diffusivity 
and L is some characteristic length. One of the difficulties in using these scalinq 
relationships is determining what length scale to use. For example, should L be 
the length, the width, the radius, the diffusion length, etc.? Since the Grashof 

number depends on the cube of the length scale, this choice can make a crucial 
difference. 


Langbein and Tiby [2] also proposed a simple criteria based on the same type of 
analysis that predicted what g-levels as a function of frequency would be required 
for different classes of proposed microgravity experiments. For an experiment in 
which heat transport is of primary interest, they estimate the relative perturbation 
due to accelerations is 
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8T 

T 


(1.4) 


§PAT 

L[co 2 + ( v / L) 2 ] 1 / 2 [co 2 + (k / L) 2 ] 172 

and for experiments where mass transport was of primary interest, the relative 
perturbation in composition is 

5C g (Ap/p) 

c l[g) 2 + ( v / L) 2 ] 1/2 [ci) 2 + (D / L) 2 j 12 ^' 5 

where g is the amplitude of the acceleration, (3 AT and Ap/p are the maximum 

relative density changes and D is the diffusion coefficient. Note that for cd-> 0 
these expressions reduce to 


5T _ g|3ATL 3 
T VK 


= Ra 


thermal 


and 


5C _ g(Ap/p)L 3 
C vD 


= Ra 


sokitao! 


which are the thermal and solutal Rayleigh numbers, respectively. 

Unfortunately, because of the difficulty in choosing the proper length scale, and 
because of the crudity of the scaling analysis, these models tended to over- 
predict the effects of residual accelerations, giving rise to unreasonable 
requirements for the acceleration environment as will be shown. 

Clearly, a better method must be found for predicting effects of steady, transient 
and periodic accelerations on various processes involving fluids. In this paper a 
powerful integral theorem is developed which, through the use of simple analysis, 
can predict the flows resulting from steady and time-varying accelerations to well 
within a factor of 2, and in many cases to within a few percent. 


1.2 DERIVATION OF THE INTEGRAL THEOREM 

The equation of motion for an incompressible fluid can be written in the form: 


dv 

p— = pvx(Vx v)-V 

at 


p+ 


pv 




+ pV 2 v + pg. 


( 1 . 6 ) 


Consider the flow around a streamline. Note that 


^vx(Vxv)*ds =|vxco»ds = 0 

since the vorticity to = Vxv is always perpendicular to v, hence Vxco will 
always be perpendicular to ds if ds is an element of length on a streamline. 
Since the pressure term (p + pv 2 /2) is a single valued scalar function, 


2 


<j> V^p + -£^- • ds = 0 

because the integrand is an exact differential. Therefore, integrating the equation 
of motion around a streamline eliminates the pressure and the inertial terms and 
yields a simplified expression [3] 

^ P lf* dS= ^V 2 v«ds + jipg.ds, Q-E.D. (i. 7 ) 

Extensive use will be made of this powerful theorem in the developing the various 
approximation formula for estimating maximum flows. 


1.3 ONE-DIMENSIONAL MODELS 

For cases in which the aspect ratio is much greater than one, it is possible to 
employ a one-dimensional model to estimate the flow field. Such models have 
the advantage of being simple enough to yield a closed form solution. 

1.3.1 Vertical Slot 

Consider, for example, the case of a fluid in a vertical differentially heated slot 
whose width is 2a and L » a. Using the Boussinesq approximation, the density 
can be written as 3 


P(x) = Po(l-P^| 

where p =(1/ p 0 )3p / 3T and x is the distance from the center line. 

Choosing a stream line that goes down at -x and up at x, integrating, and 
dividing through by p 0 , Eq. (1.7) can be written as 


^)(2L) = v^) (2 L) + aMI| (2L) 


where g is taken as positive downward. 


The boundary conditions are v(± a) = 0 (no-slip at the walls) and no net flow. Let 
^ = x/a. The steady state solution is 


3 



(1.9) 


where Gr = £ P A y w ■ w = 2a 

V 

which is the well known solution to the problem of free convection in a 

differentially heated vertical slot [4], The maximum velocity v occurs at £ = 1/ V3 
and is 


0 = 9 ft AT a 2 _ Gr v 

18V3v ” 72 V3 W ' ( 1 - 10 ) 

Note that this has the same functional form as the scaling analysis result, Eq. 
(1.1), but the maximum flow velocity is almost 2 orders of magnitude less than 
predicted by Eq. (1.1) if W is chosen as the length scale. Also note that there is 
no restriction on Gr or Pr in this model because there is no mechanism in this 
one-dimensional model for heat to be transported other than by conduction. The 
presence of ends would of course cause the flow to turn and allow for the 
transport of heat by convection, but this will be discussed in the next section. 

Also it is tacitly assumed by taking the streamline parallel to the walls that the 

Reynolds number Re = vL/ v is sufficiently small so that the flow remains 
laminar. 

The transient solution to Eq. (1.8) can be obtained by writing 
v(x,t) = v(x)+w(x,t) where w(x,t) satisfies 


9w(x,t) _ 9 2 w(x,t) 

9t 9x 2 

This has a general solution given by 

w ( x >0 = X( A n sin (b n ^)+B n cos(b^))exp(t/x n ) 

n 

where — = b 2 -^-. 

a 2 


( 1 . 12 ) 


If the initial conditions are g = 0 and v(x,t) = 0 for t < 0 and g is “switched on” at t 

= 0, the b n = n n to match the no-slip boundary conditions and the Fourier 
coefficients are evaluated to give the complete solution 


v($,t)=9^Iii 

v ’ 12v 


ft " 5 3 ) - Z ^ sin ( n 71 5) ex P0 / \ } 


(1.13) 
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The evolution of the velocity build-up is shown in Fig. 1.1. Note that the initial 
disturbance originates near the walls and the peak velocity moves toward the 
equilibrium position (£ = 1/ V3 ) as time progresses. 

A simpler approximate solution for v(t) may be obtained by integrating along the 
streamline corresponding to the steady state value of v, i. e., £, = ±1/ V3. Eq. 
(1.8) can then be approximated by v 

dv _ v gpAT 1 

dt 5 2 2 V3 O- 14 ) 

where 8 is the width of the momentum boundary layer and the negative sign was 

chosen smce the viscous drag opposes the acceleration. The solution to Ea 
(1.14) has the form 

v(t) = vH(l-e-) (,.15) 


where x = 82/v and v(<~) = ^^T 52 

v ' 2V3 v ‘ 


Since for large t, the v(t) must approach the steady state value given by 
Eq.(I.IO), S 2 can be identified as a 2 /9 and the approximate solution becomes 


v(t) = 


gpATa 2 

18V3v 


(l-e' ,/x ). 


(1.16) 


Note that the time constant x - d 2 / v - a 2 / 9 v. is very close to the time constant 
a 2 / 7t 2 v of the leading term in the Fourier expansion indicated by Eq. (1.13). It 
should be remembered that this is an approximate solution because we assumed 
the maximum velocity was along the streamline corresponding to £ = ±1/V3 , 
whereas in reality, the maximum velocity occurs at different places along the’x- 
axis as was seen in Fig. 1.1. However, this discrepancy makes little difference in 
the buildup of the transient as may be seen in Fig. 1.2. 

A closed form solution for the response of this system to periodic accelerations is 
also possible. Assume g(t) = ge w and v(x,t) = v^e** = 2A n sin(n^)e k0 '. 

Putting these expressions into Eq. (1 .8), 


ico+- 


n 2 7r 2 v^ 


A n sin(nic^) 


9 PAT t 
2 ^ 


The Fourier coefficients A n are extracted in the usual manner which yields 
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-1 


0 

x/a 


1 


Fig. 1.1 Velocity build-up as function of time for vertical slot after g is “switched 

on”. Curves correspond to t* = 0.05, 0. 1 , 0. 1 5, 0.2 where t* = t v /a 2 . The 

dashed curve corresponds to the steady state solution. Note that the velocity 
maxima is shifted slightly toward the edges at small values of t*. 
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Fig. 1 .2 Plot of v* = v(t*)/v W vs. t* = t v / a 2 . The solid line is the exact 
solution, Eq. (1.13), and the dashed line is the approximate solution, Eq. (1.16). 
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v(x) = - 


g(3ATa 2 y (- 1) n 
v n n7t(i£2+n 2 rc 2 ) 


sin(n^^) 


(1.17) 


where the dimensionless frequency £2 is co a 2 / v. For £2 = 0, this series reduces 
to the series in Eq. (1.13) with t = 0, which represents the steady state solution. 
Comparing the amplitude of the leading term in this series with the corresponding 
term in the Langbein - Tiby formula (Eq. (1 .4)), note that in this expression the 
velocity amplitude does not diminish until co is comparable with 7 t 2 v / a 2 , whereas 
in the Langbein -Tiby model, this occurs at frequencies an order of magnitude or 
more less (depending on the choice of L). In other words, the fluid responds to 
disturbances much more rapidly than predicted by the Langbein - Tiby model. 
This is compensated somewhat because the Langbein-Tiby model over-predicts 
the steady state response. 

Again a simplified approximate solution can be obtained for this system. 

Choosing the stream line at ±1/ V3 where the steady state velocity amplitude 
is maximum as the path of integration, Eq. (1.13) can be written in terms of 
maximum velocity amplitude as 


dv _ v gpAT 1 

dt 5 2 2 V3 

which, after the transient dies out, has a solution given by 

gftATS 2 1 
2V3v (ico8 2 /v+l)' 


Requiring that this reduce to the steady state solution (Eq. (1.10)) as co -> 0, we 
identify S 2 = a 2 / 9 as before, and finally 


gpATa 2 

2V3v (i£2+9) 


Comparing this to the leading term of Eq. (1.17), they are seen to numerically 
very similar. Figure 1 .3 compares these two expressions. The slight deviation at 
the higher frequencies results from the fact that the approximate calculation 

assumed the maximum velocity always occurs at £ = 1/V3, which we see from 
Fig. 1 .4 is not the case at the higher frequencies. The velocity amplitude 
obtained from the Langbein -Tiby model is also shown for comparison with L 
taken as 2 a. As discussed previously, this model seriously over-estimated the 
steady state velocity, and under-estimates the response time which can be seen 
in the figure. 
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log ( Omega ) 

Fig. 1.3 Log-log plot of the ratio of maximum velocity amplitude to steady state 
velocity vs. dimensionless frequency Q. = co a 2 / n for fluid in a vertical 
differentially heated slot. Circles represent the exact solution (Eq. (1.17)); the 
solid line the approximate solution (Eq.(1 . 1 8)), and the dashed line the prediction 
from the Langbein - Tiby model. 
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x/a 


Fig. 1 .4 Velocity amplitude scaled by steady state velocity as a function of the 
dimensionless frequency Q for the case of a vertical differentially heated slot. 
Value for Q from the top down are 0, 10, 20,50,100,1000. 
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1.3.2 Horizontal Slot 


Next consider the flow in a long horizontal slot with length L and width W = 2 a 
that is differentially heated at the end walls. A linear thermal gradient is assumed 
so that the density function may be written 

P( x ) = p 0 (1-pATx/L). 


This, of course, assumes the convective flows do not transport sufficient heat to 
disturb the linear gradient. The criteria for this will be discussed later. 

Since L»a, the viscous drag contributions from the end walls can be ignored. 

Choosing an integration path along ± y and proceeding as before, Eq (1 7) can 
be written as, 

^(2L) = v^l(2L) +g pAT( 2 y). 

Using the previous boundary conditions (no slip a the walls and no net flow), the 
second order ordinary differential equation for u has the steady state solution 

u(y) = — 6 ^ a (r|-r| 3 ) where r| = 


Again u occurs when r| = 1 /V 3 , which gives 

0 _ gftATa 3 Gr fv> 

9V3vL 72V3U, ' 


(1.19) 


This is similar to the result for the vertical slot except v/W is replaced by v/L 
Unlike the vertical slot, there is a mechanism whereby heat is transferred in a 
one-dimensional model for the horizontal slot and the above formulation is valid 
only so long as this transport does not significantly affect the isotherms. 
Otherwise the assumption that the density field is independent of the flow field is 
no longer valid. 

The transport of heat can be characterized by the thermal Peclet number 
defined as 


Pe Them,al = — (1.20) 

IV 

which may be thought of as the ratio of heat transported by convection to heat 
transferred by conduction. For the assumption that the temperature field is 
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undisturbed by the flow to be valid, the Pe Th ermai « 1. Putting Eq. (1.19) into 
this definition, we have 


p _ GrPr _ Ra Thermal 

Thermal " 72V3 = - 72 W 

where the Gr and Ra are defined in terms of the thermal gradient. Therefore Ea 
(1 . 1 9) is restricted to cases in which 

Gr p r = gp^«o(i 0 2 ). (1.21) 

(A similar restriction applies to a finite vertical slot except that Pe Therrnal would be 

given by Pe Therma | = 0 W / k . Since u is not found in the one-dimensional model 

for the vertical slot, the discussion on limits of applicability for this case must be 
deferred to the next section.) 

Examining the Langbein-Tiby model for temperature perturbations due to steady 
state convective flows, Eq. (1.4), it may be seen that, except for the numerical 

factor, the 8T/T corresponds to the thermal Peclet number if the L is taken as the 
width (or smallest dimension) of the system. 

Because of the similarity of the equations governing the flows in the horizontal 
and vertical systems, the response time of fluid in a horizontal slot will be the 
same as for the vertical slot. Therefore, Eq. (1.13), (1.15), (1.17), and (1.18) will 
apply to the case of a horizontal slot if the driving term is multiplied by W/L. 


1.4 TWO DIMENSIONAL CONFIGURATIONS 

Dressier [5] considered the flow in a vertical circular cell with a linear horizontal 
temperature distribution and obtained a transient and steady state solution in the 
that is valid for small Rayleigh numbers. Assuming v(r,0) = 0 for all r his solution 
is given by 


v e(r,t) 


gpATa 2 
32 v 


l n=1 J 0\ A n j 


exp 



( 1 . 20 ) 


where A. n are the zeros of Ji and C = r/a. 

The dominant, or most persistent term, in the series is the n=1 term which has a 
time constant given by a 2 /Xi 2 v. Since X 1= 3.8317, this time constant is 0.0681 
a 2 /v. 


This result can easily be derived using the integral condition stated in Eq.(1.7). 
For steady state the left hand side is zero. Assuming small enough Ra such that 
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the temperature field is not disturbed by the flow, the streamlines will be 
concentric circles about the origin. The viscous drag term becomes 




pd>V 2 v*ds = [i 


v + — 
v r 



27tr 


where the primes denote partial differentiation with respect to r. 
Writing g • ds = g cos0 r d8 and the density function as 

PM) = Po^-Pf^ 050 ) 

the buoyancy term becomes 

cj>p g*ds = g jjp o rcos0-|3^^-cos 2 ej de = - Tcg P ATt ' 2 


Putting these two integrals into Eq. (1.7), the differential equation for the velocity 
is 


1 dv 9 

V dt 


= v"+- 


Ve gpAT r 
r 2 4v a 


( 1 . 21 ) 


where primes represent partial differentiation with respect to r. The boundary 

conditions are “no slip” at the walls and antisymmetry about r = 0. Letting £ = r/a, 
the steady state solution can be written 


v.( r ) = ^^(C-C s ). 


which is identical to Dressler’s steady state result. The maximum steady velocity 
occurs at £ = 1/3 1/2 , which becomes 


c, _ gPATa 2 


v e = 


Gr 


48a/3v 192V3 


W 


( 1 . 22 ) 


Again the range of validity of this result is determined by the requirement that the 
thermal transport by convection be small enough so that the thermal gradient is 
approximately uniform which is required for the assumed density function to be 
valid. This requires the thermal diffusion length, 8 Therma) = k/ v 9 » a or for the 
thermal Peclet number Pe Thermal = v e a/x « 1. From the above result, we see 
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that this requires Gr v / k - Ra Them , al « 384V3 = 665 . Dressier compares his 
result to a numerical computation carried out by Robertson and Spradley [6] for a 
differentially heated square in which it is shown that Vmax is virtually independent 
of Rayleigh number up to Ra = 1000 and only falls off by 20% at Ra = 6500. For 
small Prandtl number fluids, say Pr ~ 0.01 , this model should be reasonably valid 
uptoGr ~ 10 5 -10 6 . 

An approximation for the transient can be made as before by assuming a solution 
of the form v e (t) = v e (l-e' ,/x ). Putting this into Eq. (1.21), the steady state 
portion of the solution cancels out leaving 

1 gATa 2 / = gAT 

x 32v ^ ^ 4 ^ 

From this, we find the time constant is given by 



where C is the radius of the contour chosen for the integration path. 

Fig. 1 .5 shows the velocity profiles at various times as the velocity transient 
builds up to its steady state. For steady state, £ = 1/V3 gives the maximum 

velocity and would produce % = a 2 /1 2 v = 0.083 a 2 /v which is somewhat longer 
than the dominant term in the exact solution. As can be seen in Fig. 1 . 5 , during 
the early stages of the transient, the maximum velocity occurs nearer to 

perimeter ( C ~ 1 )• Of course, taking £ = 1 would result in x = 0 since ve is zero 
around the perimeter. Selecting a C, half way between 1/V3 and 1 gives x = 
0.047 a 2 /v. Fig. 1 .6 compares the build-up of velocity predicted by these 
approximations with the exact solution. The latter estimate is seen to be a better 
approximation for the initial rise time where the higher order terms in Eq. (1 .20) 

are important, whereas some average value such as 0.065 a 2 /v gives a better 
overall fit. 


Eq. (1.21) can also be solved for periodic accelerations. Let g(t)=ge kfl ‘. After a 
transient, the velocity will take the form v(r,t)= v(r)e i<o1 . 


r^ + r*- 
dr 2 dr 


1 + 


ico r 


9 Apr 3 
2a v 


0 . 


It is convenient to express the solution as a Fourier-Bessel series, 
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Fig. 1.5 Velocity build-up in Dessler’s’ solution. Contours represent t* = 10' 3 to 
1 0-°- 3 in intervals of 1 0 - 3 . 
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t* 


Fig. 1 .6 Build-up of maximum velocity with time from Dressler’s solution (wire 
frame plot). Heavy line represents approximate solution with integration path 
taken around £ = 1/3 1/2 . Points represent the approximate solution with 
integration path taken at £ = 0.788. The dotted line represents the approximate 
solution with a time constant taken as an average of the two extreme cases. 
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1 


(1.23) 


v( r ) = _ g(Ap/p)a 2 y J,M) 

v ^ X n J 0 (2. n ) [A, 2 +i£2] 

where £2 is the dimensionless frequency, coa 2 /v. This choice of representation 
automatically enforces the no slip boundary condition at r=a. Note that in the limit 
£2 « Xn 2 , the solution reduces to 

v( r ) = _ 9(Ap/p)a 2 y 

v tt?t 3 n J 0 (X n ) 

which is a series representation of Dressler's steady state solution. On the other 
hand, if £2 » X n 2 , then 

v(r) = j g(Ap/p)a 2 y Ji(^-nC) 
n=1 X. n J 0 (X. n ) 

where it is seen that the amplitude of the velocity is inversely proportional to the 
frequency. Fig. 1 .7 shows the velocity amplitude as a function of r for various 
values of £2. The drop in amplitude does not become significant until to 

approaches the reciprocal of the time constant for the transient case, X^Wa 2 It 
is also interesting to note that for very high frequencies, the peak amplitude 
occurs closer to r = a. This is similar to the initial build up of the flow field for the 
transient case shown by Dressier. 


1.4.1 Square or Rectangular Box 

The 2-dimensional flow problem in a square or rectangular box may be 
formulated by replacing u and v in the x- and y-momentum equations (Eq. (1.6)) 

by u - ~ and v = to satisfy the continuity equation, cross differentiating, 

and subtracting to eliminate the pressure terms. The result is a single 4 th -order 

partial differential equation for the stream function y. Neglecting the inertial 
terms, this equation becomes 


+2 3V 0y_ap _ & 9p 

3x 4 ax 2 a y 2 ay 4 Mx n ay 


(1.24) 


where g* and gy are the components of the gravity vector and p(x,y) is the density 

of the fluid. The boundary conditions require that \j / and its normal derivative 
vanish at each wall to enforce the no-slip condition. 
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r/a 

Fig. 1.7 Velocity amplitudes from periodic disturbances with various 
dimensionless frequencies Q normalized by maximum steady state velocity for 
the case of a differentially heated vertical circular cavity. Value for Q from the top 
down are 0, 10, 20,50,100,1000. 
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Unfortunately this biharmonic equation has no known analytical solution. If the 
aspect ratio is very large or very small, either the x- or y-derivatives may be 
ignored to give the 1 -dimensional solutions described previously. While these 
solutions accurately describe the core flows, they cannot describe the flows near 
the ends. 

Cormack, Leal, and Imberger (CLI) carried out an extensive theoretical [7], 
numerical [8], and experimental [9] study of flows in this type of system. In their 
analytical approach, they break the model into two parts: a core region in which 
the flow is parallel and dominated by viscous effects; and end regions which turn 
the flow. For the core region, they expand the stream function, the temperature 
field, and the vorticity as a power series of the aspect ratio W/L, and find the 
solution for the stream function in the core region has the form 


xy = (C 1 +C 2 W/L + C 3 (W/L) 2 +...)(y 4 -2y 3 + y 2 ) 

where y = y/W and the C’s are functions of the Grashof number which must be 
determined by matching the core solution to the solution for the end regions. The 
problem is equivalent to the classical problem of finding the displacement of a 
hyrodstatically loaded rectangular elastic membrane with clamped edges, the 
difficulties of which are well-known. Unfortunately the solutions for the end 
regions must be carried to order 3 to get non-trivial solutions. CLI use Laplace 
transform theory techniques and obtain a solution as an expansion in Papkowich- 
Fadle eigenfunctions. The determination of these eigenfunctions must be done 
numerically, which, as it turns out, is actually more difficult than solving the flow 
equations numerically. Therefore, it is not possible to obtain a simple analytical 
expression for the stream function or the velocity profile in the end regions. They 
did, however, obtain values for the C’s in the above equation so that the velocity 
field in the core region can be written, 

U(y)= ^u ( 1_3 48x10 " Gr2pr2(w/L ) 3+ -K 2 y 3 - 3 y 2 - , -y]- 

The temperature field in the core region is given by 

T-T c x GrPr W , . ... 

TTt = L + 1335'r( 12y “S' +20y - 1 ) +0 ( Gr ^)' 

We see that the velocity will have its maximum value at y = (3 ± V3 )/ 6 which 
gives a maximum velocity of 

0 ^ (l-3.48x10^Gr 2 Pr 2 (W/L) 3 +....) 
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which is the same result obtained by Eq. (1.19) for L » W . It is interesting to 
note from this expansion that for W = L, the velocity is linear with Gr so long as 
Ra = Gr Pr « 536. This is the same order of magnitude we estimated earlier 
from the thermal Peclet number. Furthermore, from the CLI equation for 
temperature, one can deduce that the perturbation due to convective flow would 
be 


ST _ GrPr( W/L) _ gpATW 4 
T " 1440 “ 1 440 v k L 

which has the same functional form as the Langbein-Tiby model in the steady 
state limit, but numerically is diminished by more than 3 orders of magnitude. 

Bejan and Tien [10] obtained an approximation to the flow in a long horizontal 
cylindrical tube of length L and radius a with differentially heated end walls. Their 
first order solution can be written 

u (n) = s ^^(n 3 -n)s in e 


where r\ is r/a and 0 is the angular position relative to the horizontal plane. Again 
the maximum velocity will occur at 0 = ji/2 and ti = 3' 1/2 and can be written 

a _ gPATa 3 _ gpAT W 3 
12V3vL 96V3vL ' 

The flow profile in the vertical plane is quite similar to that in a horizontal slot 
except the velocity is 25% lower because of the additional drag from the 
cylindrical walls. Even so, it may be seen that the simplistic model of flow in a 
horizontal slot gives reasonably accurate estimates to the more complicated 
problem of flow in a cylindrical. 

Batchelor [1 1] considered the flow and heat transfer in a thin vertical slot 
(L»W). He expands the stream function and the temperature as a power series 
in terms of Rayleigh numbers. For the first order term, valid for small Ra and L - 
W , he uses an approximate stream function which he attributes to Grashof (no 
reference given) which can be written, 

(x } gP(AT/L)W 4 2 ( (x/L) 2 (1- x / L) 2 (y / W) 2 (1 - y / W) 2 
v 3 1 (1+W 4 /L 4 ) 

where L is taken along the horizontal or x-axis. 

For the case of a square, the velocity may be found from 
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By _ 2 
By - 3 


Gr^(x/L) 2 (1-x/L) 2 [2(y/W) 3 -3(y/W) 2 +y/w]. 


(1.26) 


The maximum u will occur at x = L/2 and y = (3± V3 )W/6 which becomes 
~ _ Gr v 

U “ 144V3 W ' < 1 - 27 ) 

This model, even though it is not an exact solution to the biharmonic equation, 
produces streamlines and velocities that agree reasonably well with numerical 
computations for small Rayleigh numbers as can be seen in Fig. 1.8. 

In order to apply of the integral theorem to flows in square or rectangular cavities, 
one is tempted to use the functional form of Eq. (1 .25) to evaluate the derivatives 
in the viscous term. However, this leads to an inconsistency as will be shown. 
For example, if one chooses the perimeter of a square cavity as the streamline 
for the path of integration, the contribution from one side is 



,2,,N 


a z u 

By 2 


/ y=0 


dx . 


It is convenient to combine Eq. (1 .25) and (1 .26) and write 


u = -^u(x/L) 2 (1-x/L) 2 [2(y/W) 3 -3(y/W) 2 + y/w] 


(1.28) 


The contribution from the B 2 u/Bx 2 term vanishes since u is identically zero 
everywhere along the path. The contribution from the second term yields 


r B 2 u 288 1 -6 ^ 

] By 2 ° X " V3 30 v W 2 J 


Lu. 


For L - W, the other 3 walls will give the same contribution since u and v are 
equal in magnitude. The driving force is given by 

<fp 0 (1-|3ATx/L)g.ds = p 0 f3ATL. (1.29) 

Equating the terms gives a maximum velocity 

Q= 15 gftATW 2 _ 15 Gr v 
8 144-V3 v 8 144V3 W' 

This is inconsistent with the maximum velocity obtained from the velocity function 
Eq. (1.27), which is indicative of the fact that, even though the velocity function in 
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0.0003 



Fig. 1.8 Comparison of stream function computed from Eq. (1.25) (dotted curve) 
with numerical computation for L = W. The quantity g p AT/v was taken as 
0.1915. 
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Eq. (1.26) gives reasonable results, it does not exactly satisfy the momentum 
equation. Or for that matter, neither does the stream function given in Eq. (1.25) 
exactly satisfy the biharmonic equation as can be verified by differentiation. This 
can be seen in Fig. 1 .9 and 1 ,9A which compares the stream function computed 
numerically with Eq. (1 .25) in the vicinity of y = 0. Even though the agreement is 

close, the third derivative of 'F computed numerically is slightly more than twice 
that obtained from differentiating Eq.(1.25), which would explain the over- 
prediction of the velocity by almost a factor of 2. As it turns out, the damping 
term computed by differentiating Eq. (1 .25) at x = L/2 is a very good 
approximation for the average damping along the integration path. Therefore, 



Each of the other 3 walls will give the same contribution, hence the maximum 
velocity becomes 


fl= g|3ATW 2 Gr v 
144V3v 144V3W 

which agrees with Eq. (1.27). 

As was pointed out by Batchelor, extension of this model to rectangular 
configurations quickly breaks down as the aspect ratio departs from unity 
because the vertical flow field is unrealistically spread over half the length of the 
channel. This is seen in Fig. 1.10-1.12 which compares the x-dependence of the 
stream function along the center line computed from Eq. (1.25) with the 
numerical results. 

Experience from the numerical modeling suggests a model for rectangular 
configurations in which the end flows turn within the region L/2 for W < L < 2 W, 
and within W for L > 2 W, assuming that L is taken along the horizontal or x-axis. 
Using this model, Eq. (1.25) can be differentiated to obtain the maximum 
horizontal and vertical velocities for L < 2 W which yields 


1 1 — 


Ll — 

^ L 

from which 



Gr v 

u = 

1 

72 V3 L 


_ Gr v fy/W-3(y/W) 2 +2(y/W) 3 


12 L 


I 


1 + W 4 / L 4 


2 L> W 


(1.29) 


(1.30) 


and 
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0.0003 



curve 


Fig . 1 .9 Shape of the stream function computed numerically compared to Eq. 
(1.25) (dashed line) normalized to give the same peak value. The numerical 
values (circles) were fit with an 8th degree polynomial (solid line) which was 
differentiated at x = 0 to obtain 'f'”’(0) = 0.10972, slightly more than twice the 
value 0.0496 from Eq. (1.25) even after it was normalized to give the same peak 
value. 
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Fig. 1.10 Comparison of stream function computed from Batchelor’s approximate 
model (dashed curve) with numerical computation for L = 2 W. 
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x/W 

Fig. 1.1 1 Comparison of stream function computed from Batchelor’s approximate 
model (dotted curve) and from Eq. (1.33) (dashed curve) with numerical 
computation (solid line) for L = 5W. 
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Fig. 1.12 Comparison of stream function computed from Batchelor’s approximate 
model (dotted curve) and from Eq. (1.33) (dashed curve) with numerical 
computation (solid line) for L = 10W. 
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v = 2t) _ G L vWfx/L-3(x/L) 2 +2(x/L) 3 l 

^xJy=w/ 2 12 L L j 1 + W 4 /L 4 J ^* 31 ) 

from which 


72V3 L L v 1 + W 4 /L 4 J L< ^ VV - 0-32) 

Since the magnitude of the stream function peak is predicted with reasonable 
accuracy by Eq (1.25) as was seen in Fig. 1.10-1.12, Eq. (1.29) will give a qood 
i a A? P I° Ximatl ? n for the maximum horizontal velocity for all L > W/2. In fact for L » 
Z' Fn S n ,n m V ? dU f e ? o^'mensional result for the horizontal slot given 

ma^i mirtiSr 3 ^ 1 L < W - (130) gives a good approximation for the 
1 4/wT U ?L ert Ca v '® ,oclty as can be seen by multiplying the top and bottom by 
L / W • This gives the same form as Eq. ( 1 .29 ) except L and W are 

imerch a nged This reduces to the one-dimensional vertical slot Eq. (1.10) in the 
limit L « W if the length scale in the Gr is taken as the smallest dimension - in 

Inis CciS0 L. 

To obtain a reasonable estimate for v for L > 2 W, we introduce a stream 
function for the end region that peaks at W using a form similar to Eq. (1.25), 


u/r v i/ox- 9P(AT/L)W 4 1 f (x/2W) 2 (1-x/ 2W1 2 
v 24 [ (1 + W 4 /L 4 ) 

This can be differentiated to give 


(1.33). 


Gr \f 1 A 
1 44 V3 L v 1 + W 4 / L 4 J ’ 


L > 2 W 


(1.34) 


Similarly to obtain u for W > 2 L, we introduce a stream function that peaks at 


y ( 1/ 2 j y) _ gP(AT/L)L 4 1 I (y / 2L) 2 (1 — y/ 2 L) 2 
v 24 [ (1 + L 4 /W 4 ) 

which can be differentiated to give 


(1.35) 


Gr vW ( 1 


Gr* v 


144V3L L v 1 + W 4 / L 4 J 144V3 L(iTL 4 / W 4 J ’ W>2L 0-36) 


where Gr* = g p AT L 3 / v 2 . 


29 



To obtain an approximation for the time constant for response to transient 
accelerations, the integral theorem may be used. Choose horizontal streamline 
along + a’ relative to the center line and vertical streamlines along ± b’ and use 
the average values for viscous damping terms as described previously. Thus 
Eg. (1.7) may be written 



= 1^1-3 AT^jg.ds 


To enforce the continuity equation, derivatives of the stream function will be used 
in lieu of velocity functions. Inserting these and performing the indicated 
integrations, we obtain 


WX W'/d'fM 
k ay/ + L' \3x/ b , + v 


a 3 '? 


y=a' 


+ v\r /at\ 

+ u \ ax / h . 

' / x=b 


g(3AT W' 
2 L' 


For L = W, the derivatives of V may be obtained from Eq. (1.25). Assume *F = 

'F (1 - e ^ T ) where *F is the steady state solution to the above equation. Inserting 
this function, the time constant may be obtained as 

t = W 2 [ (2y 3 -3y 2 + y)+(W/L)(W7L , )(2x 3 -3 x 2 + x)| 

X v { (l2y-6) + (W/L) 3 (W7L')(l2x-6) J (1 ‘ 38) 

where W’ = 2 a’, L’ = 2 b’, and x and y are the values of x/L and y/W 

respectively for which the velocity is maximum. For L = W it is reasonable to 
assume that L’ and W’ stand in the same ratio as L and W. For steady state, the 

maximum velocity is obtained for x = y = (3- V3)/6. Putting this into the 
above, we obtain 


x 



1+(W/L) 2 
36[l+(W/L) 4 ] ' 


(1.39) 


For W = L, this reduces to x = W 2 /36 v = a 2 / 9 v, the same as for the one- 
dimensional case. If L » W, or W » L, this also reduces to the one-dimensional 
case with the shortest dimension as the length scale. 
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The maximum horizontal and vertical velocities computed from the above model 
are compared with numerical results for a variety of different shapes in Table 1. 
Also Fig. 1.10-1.12 compares the model stream function along the center line 
with numerical results. It can be seen that the model is guite accurate for square 
or large L/W ratios but is somewhat less accurate for the intermediate cases. 
This is to be expected considering the crudity of the approximation that 
describes the turning points in the flow. However, even in the worst cases the 
model is off by only 20%. 


1.5 FLOWS AT LARGE RAYLEIGH NUMBERS 

Thus far we have considered cases where the Ra was small enough so that the 
isotherms were unperturbed by the flow. Let us now expand the approximation 
method to the case of large Rayleigh numbers where the heat transfer is 
dominated by the convective flow in the thermal boundary layer. For the time 
being we ignore the horizontal and concentrate on the vertical flow. The heat 
transfer equation must also be solved with the flow equation. For the case at 
hand the terms of interest are 


v( x )£L - K 
[ By dx 2 


Approximate ~ and where 5j is the thermal length or the 

distance over which the temperature changes significantly due to convective flow. 
If this length is long compared to L/2, heat transport is dominated by conduction; 
if the thermal length is short compared to L/2, convective heat transfer 

dominates. Putting these back into the heat transfer equation, the &r may be 
estimated as 


1 __ v 

5* “ kW ‘ C 1 - 40 ) 

The viscous damping integral may be written 


j>V 2 v«ds 



9 3 'f / (x,1/2) 

dx 3 


x=0 


dy+2 


f/ a 3 y(i/2, y ) 

J o\ ay 3 



dx 


(1.41) 


We assume the flow field is contained in some thermal boundary layer b, which is 
proportional to 5j, and can be expressed as 


v(x) 



(2X-3 X 2 



X = x/b . 
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Table 1.1. 



W 

(cm) 

g P AT /v 
(cm- 1 sec- 1 ) 

Ucomp 

(cm\sec) 

v comp 

(cm/sec) 

U e st 

(cm/sec) 

Vest 

(cm/sec) 

i 

1 

0.1915 

7.94x1 0- 4 

7.94x10-4 

7.68x10-4 

7.68x10-4 

2 

1 

0.1915 

7.27x10-4 

4.50x10-4 

7.22x10-4 

3.61x10-4 

5 

1 

0.1915 

3.26X10- 4 

1.64x10-4 

3.00x10-4 

1.50x10-4 

10 

1 

0.1915 

1.56x10-4 

9.68x10-5 

1.54x10-4 

7.58x10-5 

5 

1 5 

0.0050 

4.96x10-4 

4.96x10-4 

5.01x10-4 

5.01x10-4 

1 

5 

0.0050 

2.50x10-5 

4.05x10-5 

1.96x10-5 

3.91x10-5 

10 

10 

2.83x10-5 

1.13x10-5 

1.13x10-5 

1.14x10-5 

1.14x10-5 

2 

10 

2.83x10-5 

5.19X10- 7 

9.79X10- 7 

4.43x10-7 

8.86x10-7 

10 

2 

2.83x1 0- 4 

1.81X10- 6 

1.04x10-5 

1.77x10-5 

8.86x10-7 


Comparison of numerical computations of maximum horizontal and vertical flow 
velocities with approximates for small Rayleigh Numbers. L is taken to be 
perpendicular to the g-vector. 


32 


























































The stream function as a function of x may be obtained by integrating this 
expression, 


^(x.1/2) = 4'F( X 2 -x 3 + % 4 /4) where * = ^-vb. 

' 8 

Putting this into Eq. (1 .40), we get 


1 

6? 


8*F 

3V3 bxW 


(1.42) 


We relate b to 5j by b = y 6j where y is a constant of proportionality to be 
determined. Putting this in the above expression we get 


8V 


b 3V3 KWy 2 


(1.43) 


Differentiating the stream function Eq. (1.42), we obtain 
d 3x F(x)"| _ 24 V 

a * 3 


(1.44) 


which, with the help of Eq. (1.43), can now be written as, 


9 3 ^(x) 

3x 3 



24 -8 3 V* 

(3V3)Y(^WT 


(1.45) 


Ignoring the other terms for the time being and assuming this term represents the 
average along the path of the integration as before, the maximum stream 
function can be obtained by integrating Eq. (1.41) and equating it to the driving 
term which is still given by Eq. (1 .29) 


24 -8 3 y 4 _ gpAT 

(3a/3)V (kW) 3 2v 

from which 


8 

(1.46) 


^V3 V " 


, 3/2 rgpATW 3 K 3 


, 1/4 


= 0.275 y 3/2 (GrPr) V4 K 




The maximum velocity may be obtained from the relation 
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1/2 


(1.47) 


A _ &± 1 _ f 8 Y * 2 
3V3b UV3j k Wy 2 




v 

W 


A series of numerical computations were carried out with large Ra so that the 
flow would be dominated by the thermal boundary layer at the walls and the 
maximum stream function and maximum vertical velocity was correlated with the 
above results (see Table 1 .2). The choice of y is determined by attempting to 
obtain a consistent model for both the maximum velocity and the maximum 
stream function. The velocity data would argue for a value of -1.3 whereas the 
stream function data would argue for -1.7. Therefore, a compromise choice of 

1 .5 is taken. The fact that different values of yare required to correlate the 
different quantities reflects the inability to describe the flow field by a simple one- 
dimensional expression representing some average flow. However, the fact that 
the range of y needed to correlate the velocity and stream function data is 
relatively narrow does indicate that despite the crudity of the model, it does seem 
to have useful predictive capabilities. 


Now an interpolation equation for the transition between low and high Ra can be 
developed. Assume some sort of effective boundary layer in which the flows 
along the vertical faces are contained. This effective boundary layer should 
approach b at large Rayleigh numbers when flows are dominated by the thermal 
boundary layer and should approach L/2 (or W, whichever is smaller) for smaller 
values of Ra where the isotherms are not carried by the flow. There are several 
ways of representing a function with these properties. For example, one could 
define 


1 


5 


eft 



L <2 W , 


(1.48) 


or one could just as easily define 


a - eKe 1 ; L<2W ' 


J e« 


(1-49) 


or something in between. 


We first choose the form of Eq. (1 .48). Now Eq. (1 .44) becomes 


Using Eq. (1.43) for 1/b, this becomes 


34 



Table 1 .2 


Table 1.2. 

K 

(cm 2 /sec) 

V 

(cm 2 /sec) 

Pr 

Gr 

Ra 

T comp 

(cm 2 /sec) 

v comp 

(cm/sec) 

Yi 

Y2 

0.001 

0.010 

10 

625000 

6250000 

0.031 

0.117 

1.70 

1.31 

0.001 

0.010 

10 

62500 

625000 

0.018 

0.0378 

1.73 

1.34 

0.001 

0.010 

10 

6250 

62500 

0.010 

0.0116 

1.70 

1.30 

0.001 

0.010 

10 

625 

6250 

0.004 

0.0030 

1.41 

1.05 

0.010 

0.010 

1 

625000 

625000 

0.162 

0.349 

1.64 

1.23 

0.100 

0.010 

0.1 

625000 

62500 

0.718 

0.850 

1.40 

0.95 

1.000 

0.010 

0.01 

625000 

6250 

3.530 

2.20 

1.28 

0.78 


Computations at large Ra used to determine the factor y . Maximum vertical 
velocities and stream functions were computed numerically for a 5 cm x 5 cm 
chamber with a Gr = 6.25 x 10 5 with different values of k and v. The yi is the 
parameter needed to reconcile Eq. (1.46) with the numerically computed 
maximum stream function and ^ is that needed to reconcile Eq. (1.47) with the 
numerically computed v. 
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3x J,.0 S’„ 3V3 kW y ! lJ 

The y-derivative can be expressed as 


L < 2 W 


a 3 v( y) 
a y 3 




Taking these two derivatives as the average terms in Eq. (1.41) and equating the 

viscous term to the driving term using y = 3/2, a 4th order polynomial for *F is 
obtained, i.e., 


'T 


/ 8~2 * | W]' f L gpATW 3 
y9■3^f3 k L J W 2 3 48 v 


or 


*F 


16 W W 
27^3 k + L 


, xp L _ gpATW 3 
W 384 v 


Had we chosen the for of Eq. (1 .49), this would be 




16 'F 
27V3 






+ y_L - Grv 
W 384 


(1.50) 


(1.51) 


These relations hold for L < 2 W. For L > 2 W, replace the W/L term by 1/2. 
Note that for W = L and 'F « k , these reduces to 


y = gftAT(W/2) 3 
96 v 


which is the maximum value for Eq. (1.25) with W = L. For *F » k, Eq.(1.46) is 
recovered. 

Having found the maximum stream function, the maximum horizontal and vertical 
velocities may be determined from, 


- _ 8»F 1 
3V3 a ’ 


8*F 1 _ 84* f 32 T 2 > 
3V3 6 e(f " 3V3[27V3 kW + L / 


(1.53) 
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Again, for L > W/2, substitute 1 for 2/L. 

This model was tested against a series of numerical computations for W = L over 
a wide range of Ra and k and the results are shown in Fig. 1.13-1.15. According 
to Eq. (1.50) or (1.51), the dimensionless quantity obtained by dividing the stream 
function by k should be a function only of Gr Pr = Ra and the aspect ratio. 
Similarly, the horizontal and vertical velocities can be expressed in terms of 

thermal Peclet numbers (u or v times W/k) which again should be functions only 
of Gr Pr and the aspect ratio. 

As may be seen in Fig. 1.16, the stream function determined from Eq. (1.51) 
seems to give better agreement with the numerical results for Pr > 1 , whereas 
Eq. (1 .50) seems to give better agreement for Pr <1 . A similar trend is seen in 
the comparison of vertical velocities in Fig. 1.14, although the model slightly over- 
predicts the velocity for the small Pr cases at large Ra. Actually, there was a 
convergence problem with the numerical computations for Pr = 0.01 at high Ra, 
hence the velocity data for these points may not be accurate. Furthermore, the 
Reynolds numbers for these flows are 0(1 0 3 ) which means that they may be 
nearing a bifrucation point where the flow is no longer laminar. The weakest 
feature of the model appears to be the prediction of horizontal velocity for Pr > 1 
at high Ra as may be seen in Fig. 1.15. This comes from the assumption that the 
y-dependence of the horizontal was still give by Eq. (1 .29) for large Gr Pr, which 

is probably reasonably valid for L » W, but not for L = W. 


1.6 EFFECTS OF AXIAL DENSITY GRADIENTS 

Thus far, it has been assumed that the g-vector was always perpendicular to the 
density gradient. However, most systems of interest will have a thermal of solutal 
density gradient along the axis of symmetry and the g-vector can have both axial 
and transverse components. Therefore, it is necessary to investigate the 
possible effects of an axial stabilizing or destabilizing density gradient on the 
flows produced by a transverse acceleration. 

Consider a two-dimensional chamber with length L and width W oriented nearly 
along the g-vector but with a slight tilt. Take the x-axis along the length and the 
y-axis in the transverse direction so that the g-vector can be resolved into 
components g x and g y with g x » g y . Assume a linear thermal gradient is 
imposed along the x-axis. 

As stated previously, neglecting inertial effects, the steady state flow can be 
described by Eq. (1 .24) 

V 4 'F = MZL_M^I 

v dx v 3y 
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log (Ra) 

Fig 1.13. Maximum stream function normalized by the thermal diffusivity vs. Ra 
for different Prandtl numbers. Symbols represent Pr = .01 (diamonds), Pr = 0.1 
(triangle), Pr =1.0 (square), Pr = 10 (circles). The solid and dashed line 
represent two different interpolation schemes in the approximate model (see 
text). 
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log (Ra) 


Fig 1.14. Vertical Thermal Peclet number vs. Ra for different Prandtl numbers. 
Symbols represent Pr = .01 (diamonds), Pr = 0.1 (triangle), Pr =1.0 (square), Pr = 
10 (circles). The solid and dashed line represent two different interpolation 
schemes in the approximate model (see text). 
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Fig 1.15. Horizontal thermal Peclet number vs. Ra for different Prandtl numbers. 
Symbols represent Pr = .01 (diamonds), Pr = 10 (circles). The solid and dashed 
line represent two different interpolation schemes in the approximate model (see 
text). 
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log (Pr) 


Fig. 1.16 Plot of the Reynolds number squared vs. Pr for constant Gr = 6.25 x 
10 5 . The solid line represents the scaling law Re 2 ~ Gr/Pr, the dashed line 
represents the predictions based on Eq. (1.53), and the dotted line represents 
scaling based on Goldstein’s approximate solution for a free heated vertical wall 

in which Re 2 ~ Gr/(Pr +20/21). Both scaling laws were adjusted to converge at 
Pr =10. 
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This formulation assumes that p = (1/p)3p/5T<0 and g x and g y are taken as 
positive in the x- and y-directions respectively (or that P = - (1/p)dp/3T>0 and 
the accelerations are taken as positive along the positive coordinate directions). 

The heat flow is described by 


dT dT f d 2 T a 2 T^| 
u ax + v a y K (ax 2 + a y 2 / 

Since the core flow is essentially one-dimensional, Eq. (1.24) reduces to 


(1.54) 


V 4 V _a 3 u _g y paT g x paT 
ay 4 ay 3 v ax V ay 

Assume the temperature field can be described 


by 


(1.55) 


T(x,y) = Y^x+f(y) 
which allows Eq. (1.54) to be reduced to 


AT 

u t =k 


f a 2 ?" 1 


a y 2 


(1.56) 


(1.57) 


Introduce a dimensionless length r\ = y/a where a = W/2. Eq. (1 .55) and (1 .57) 
can be written 


U '" (T1 ) = QxPT'a 2 | g y p(AT/L)a 3 

V V 


(1-58) 


and 


u(AT/L)a 2 = kT" (1.59) 

where the prime denotes differentiation with respect to t]. Differentiating Eq. 
(1.58 again and using Eq. (1.59) to eliminate f" gives 

u^Cti) = u(ri) =-Ra u(q). (1.60) 

The solution to this equation is 

u(T|) = C 1 sinh(KTi) + C 2 cosh(KTi) + C 3 Sin(KTi)-i-C 4 cos(K'n) (1.70) 
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where K = (-Ra) 1/4 . Conservation of flow requires J u(r|)dT| = 0 which eliminates 

-1 

the even terms; hence, C2 = C4 = 0. The “no slip” boundary condition at the walls 
requires 


C 1 sinh(K) + C 3 sin(K) = 0. (1-70) 

Putting the velocity function u(q) from Eq. (1 .70) back into the heat transfer 
equation, Eq. (1.59), yields 


(C 1 sinh(KTi) + C 3 sin(Kri))^I^ = K f". 
Integrating twice gives 


(1.72) 


T(ti) = 


_ ATaVQ, . 


C, _. 


Lk ^ sinh ( KT l)~j^ sin ( KT l) + C 5 T lJ- 

Putting the first integral of Eq. (1 .72) back into Eq. (1 .58), we obtain 

g y P AT a 3 g x p AT a V C, u/ „ x C, _ 

( ^ ) = ~ l v ' k L [~ k ' cosh ( KT l)-jf cos ( KT l) +C 5 

= C, K 3 cosh(K r() - C 3 K 3 cos(K q) 
which allows C5 to be identified as 


(1.73) 


(1.74) 


C 5 


9y K 

9x a' 


This allows the temperature to be written as 


(1.75) 


T(t|) = 


AT a 2 
Lk 


§-sinh(Kri)-^f-sin(KTi)- 


•^1 

g x a ) 


(1.76) 


The remaining coefficient must be determined from the thermal boundary 

condition at the walls. For insulating walls, f'(±1) = 0. Differentiating Eq. (1.76), 
we obtain 


C, cosh(K) - C 3 cos(K) = K 

g x a 


(1.77) 
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which together with Eq. (1.71) allows the determination of Ci and C 3 . The final 
result can be written 


u(ti) = 


K 9y_ k ( -sin(K)sinh(Kr|) + sinh(K)sin(KT|) ^ 
g x al^ sin(K)cosh(K) + cos(K)sinh(K) } 


and 


T(ti) = 



9y _ af sin(K)sinh(Kr|) + sinh(K)sin(Kr|) 
g x K^ sin(K)cosh(K) + cos(K)sinh(K) 



(1.78) 


(1.79) 


For a stabilizing gradient, the Ra is > 0 , which means that K is complex, i.e., 


K = (-1) ,/ 4 Ra 1/4 = 


2 2 


Ra = y(1 + i). 


(1.80) 


Putting this back into Eq. (1.78, we get 


U(n) = 2Y^- 
9x a 


cos(y)sinh(Y)sin(yr|)cosh(YTi)-sin(Y)cosh(y)cos(yr|)sinh(Y-n) 

sin(y)cos(y)+sinh(y)cosh(y) 


(1.79) 

This may be put into a more recognizable form by expanding about yand then 
replacing y with V2Ra 1/4 / 2 to get 


U(rj) = 9yP( AT/L ) a 
v 


krf - n) - (r | 7 - 7 Ti 5 + 35 if - 29 ti ) + 0(Ra 3 )] (1.80) 


The first term in the series represents the core flow in a horizontal slot and the 
second term represents the first order correction from the stabilizing gradient 
characterized by Ra. However, note that the correction is minuscule for small 
values of Ra for which this expansion is valid. 


Similarly, the temperature perturbation may be written 


T(ti) = sin(y)cosh(y)sin(yr|)cosh(yT|) + cos(y)sinh(y)cos(yT|)sinh(yr|] 

sin(y ) cos(y ) + sinh(y) cosh(ri) 


l gj 


Expanding about y and replacing y with V 2 Ra V4 /2 as before 


(1.81) 


-yn 
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f(t|) = AT-g y - t3(AT/)a4 
L 1 080 k v 


(1.82) 


(9r] 5 -30ri 3 + 45Ti) + O(Ra 2 ). 


(Note: This is identical to the result of Cormack, Leal, and Imberger discussed in 
Section 1 .4.1 as may be seen by setting y = 0 and W = 2 a in their expression 
and comparing it to the above expression with r| = 0.) 

Since the coefficient of Eq. (1 .78) can also be written as 

9y k = g y P(AT/L)a 3 = 9V3 u 0 
g x a Ra(x)v Ra(x) 

where u 0 is the maximum velocity that would be attained in the absence of the 

effect of the x-component of acceleration which can be obtained from Eq. (1 .80), 
it can be seen that the ratio 


u(tQ _ 9V3K r -sin(K)sinh(Kr|) + sinh(K)sin(K'n) ^ 
u 0 Ra(x)[ sin(K) cosh(K) + cos(K) sinh(K) } 


(1.83) 


is a function of Ra(x) only. Fig. 1.17 shows this velocity ratio for various 
stabilizing and destabilizing values of the Ra. 

Similarly, since the maximum temperature perturbation from the y-acceleration 

driven flow form Eq. (1.82) is seen to be f 0 = Ra M ATW , the temperature 

45 L 

perturbation can be normalized to 


T(r|) _( 45 1 T sin(K)sinh(Kr|) + sinh(K)sin(KTi) 

f 0 V Ra ( x ) k( sin(K)cosh(K) + cos(K)sinh(K) ^ 


(1.84) 


which is also a function of Ra(x) only. Fig. 1.18 shows how this ratio varies with 
various stabilizing and destabilizing values of Ra. 

Notice that the velocity and temperature perturbation decreases slowly with 
increasing Ra > 0, but increases rapidly as Ra < 0. This is seen more vividly in 
Fig. 1.19. In fact there is a singularity at K = 2.365 when the denominator of Eq. 
(1 .79) and (1.81 ) goes to zero. This corresponds to Ra(x) = -31 .285, which 
presumably is the critical Ra for this configuration. Inspection of the 

corresponding equations for Ra > 0 where K = y (1+i) shows no singularities. 
The normalized effect of an axial acceleration can be represented empirically by 
IRacl/(IRa c l+Ra(x)) as may be seen in Fig. 1.19 From this, it is apparent that the 
stabilizing Ra(x) must be substantially greater than IRa c l in order to be effective. 
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Fig. 1.17 Core velocity for different axial accelerations normalized by maximum 
core velocity with no axial acceleration (Ra=0). Solid line corresponds to Ra = 0; 
dot-dashed line to Ra = -10; dash-dash-dot line to Ra = -20; dashed line to Ra = 

10; dotted line to Ra = 100. Here Ra = 9*fl( AT/r| -) a4 

KV 
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Fig. 1.18 Temperature perturbation for different axial acceleration normalized by 
the maximum perturbation for no axial acceleration (Ra=0) Here 
Rn . g x P(AT/L)a 4 

KV 
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Ra (x) 

Fig. 1.19 Normalized velocity taken at t \ = -1/3 1/2 (solid line) and normalized 
temperature perturbation (taken at t \ = 1) as a function of Ra(x). The dotted line 
is the empirical function IRa d /(Ra+IRa C |). Here Ra = g x P(AT/L)a 4 / kv and 
Rac =-31 .285. In terms of the more conventional definition 
Ra = g x P(AT/L)W 4 /Kv, Rac = -500.56. 
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Therefore, it does not appear to be practical to attempt to use an axial 
acceleration to stabilize microgravity experiments involving dilute systems with 
low Pr fluids against unwanted low-level transverse accelerations because the 
Gr(x) required to achieve a large enough Ra(x) to be effective would no longer 
qualify as a microgravity experiment. There is still considerable benefit to be 
derived by attempting to orienting the furnace axis along the residual acceleration 
vector, if for no reason other than to reduce the more critical quasi-steady 
transverse accelerations. 

Non -dilute systems in which the solutal gradient is much larger than the thermal 
gradient can also be analyzed with this model by substituting compositional 

change AC for AT, Ap/p for (3AT, the diffusion coefficient D for k, and the Schmidt 
number Sc for the Prandtl number Pr. Thus the thermal Ra, which Gr Pr 
becomes the solutal Ra which is Gr Sc. Since Sc are generally much larger than 
Pr for systems of interest, considerably more benefit might be realized by 
applying a small stabilizing axial acceleration to reduce the effects form 
uncontrollable transverse accelerations. 

The analysis can also be applied to a vertical Bridgman system with a slight tilt in 
unit gravity. Here the Ra - 0(10 4 -10 5 ) and the result is shown in Fig. 1.20 and 
Fig. 1 .21 . As the stabilizing field is increased, the residual flow is more or less 
confined to the vicinity of the walls and the temperature perturbation, when added 
to the linear unperturbed gradient, forces the resulting isotherms (or iso- 
concentrates) to be essentially horizontal with respect to the resultant g-vector. 


1.7 EFFECTS OF MAGNETIC FIELDS 

Unwanted flows in solidification process are often reduced to tolerable levels on 
earth by the application of strong magnetic fields, either axially of transversely to 
the melt. Similar benefits can be achieved in microgravity experiments as will be 
shown. To analyze these effects, it is important to have a more accurate 
description of the flow in the end regions than could be obtained by the crude 
assumption suggested earlier in Section 1.41. For configurations in which L » 
W (here L is taken as the longest dimension), we have shown that for pure 
viscous damping the flow in the core region is given by a simple cubic function of 

the form r| 3 - q where r| is the distance from the centerline normalized by a, the 
half width. If we assume this relation holds over the entire length, we can 
approximate the stream function by 

'F(^ri) = 0(^)(r| 4 — 2q 2 + l) (1.85) 


where £ is the distance along longest dimension normalized by the half-width, a. 

Inserting this into Eq. (1.24) and performing the necessary differentiations along 

the center line (r) = 0), one obtains an ordinary fourth-order differential equation 
of the form 
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-1 0 1 

y/a 


Fig. 1.20 Normalized velocities For Ra = 10 4 (solid line), 10 5 (dashed line), and 
10 6 (dotted line). Note that as the Ra increases, the residual flow tends to 
essentially creep along the walls. 
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Fig. 1 .21 Normalized temperature plot for large Ra. Note the perturbation is 
nearly a straight line except at the walls. When this perturbation is added to the 
added to the linear unperturbed gradient, it forces the isotherms to be essentially 
perpendicular to the resultant g-field. 
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<t>"" - b O" + c <I> = G 

where b and c are constants. The G is identified as 


d P 


(1.86) 


G t* uu 

= — -9y^: + g 


ax 3x a y 


Using Laplace transforms with initial conditions 0(0) = <D'(0) = 0, the solution 
may be written 

= r^ 2 + r a g [C 1 (cosh(R 1 ^)-cosh(R 2 ^))+ 

+ C 2 fAsinh(R^)-^sinh(R 2 ^)l+fJ T cosh(R^)--Lcosh(R 2 ^)] 

V n i n 2 ) ) 


*'&)= 


Q 

= R~ 2 -R 2 2 [ Ci ( Ri sinh ( R r ^)- R 2 sinh(R 2 ^))+C 2 (cosh(R 1 ^)-cosh(R 2 ^))- 
+ f^-sinh(R^)--L S inh(R 2 ^)l 

V n i M 2 )_ 


where R1 and R2 are the roots of the indicial equation 


R - b R + c = 0 (1.9( 

given by 

R 1 = ^-V b2 + 2 c 2 + ^-V b2_2c2 and R 2 = V ^ 2 + 2 c 2 -j-^/b 2 -2c 2 

(1.9 

If 2 c 2 > b 2 , which will generally be the case, the roots will be complex. 

The coefficients Cl and C2 are determined by the remaining boundary 
conditions, i.e. 0(L/a) = <I>'(L/a) = 0. 

The coefficients may be evaluated from 


C — ^ 22 ~l~ ^2 ^12 j Q _ B 2 A„ + B 1 A 21 

' A A A A 2 A a x a 


^11^22 ^21^12 


An ^22 ^21 ^12 


where 
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A„ = cosh(R 1 L/a)-cosh(R 2 L/a) 

A 12 = ^-sinh(R 1 L/a)-^-sinh(R 2 L/a) 

Ri R 2 

A 21 = R 1 sinh(R 1 L/a)-R 2 sinh(R 2 L/a) 

A 22 = cosh(R 1 L/a)-cosh(R 2 L/a) 

B i = + cosh ( R i L/a )-^rCosh(R 2 L/a) 

B 2 = ^-sinh(R 1 L/a)-^-sinh(R 2 L/a) 

Ri R 2 

If <!>(£,) is more or less constant away from the end regions, the solution may be 
simplified considerably by using the semi-infinite approximation; i.e., <I> and O' 

must remain finite as x -> «. To enforce this condition, the hyperbolic functions 
with their complex arguments were expanded into real and imaginary parts and 
the coefficients Cl and C2 were chosen so that the positive exponential terms 
vanished. Writing the complex roots = a + i p and R 2 = a - i p so that a = (Ri 
+ R 2 )/2 and j3 = i (Ri - R 2 )/2, the resulting expression is obtained; 


= G L 


^sin(p^)+cos(p^)je- a5 


a +2 a P +p 


G 


1 -^sin(P^)+GOs(p^) 




(1.92) 


In cases involving large aspect ratios and where the semi-infinite approximation 
is not appropriate, it is more convenient to take the origin at the center of the slot 
and require that <I>(L/2) = <I>'(U2) = <t>(-L/2) = 0'(-L/2) = 0. Now the solution can 
be written 


<&(£) = 


G 


R, 2 R 2 2 R 1 2 -R 2 2 


f A 

C 1 + C 2 R 2 +-1 2 
M 1 


cosh(R 1 £)-[c i +C 2 R 2 2 +— ! y cosh(R 2 £) 
V R 2 J 


(1.93) 


*'(5) = 5 T^|fc,R,+C,R, 3 +J-)slnh(R,5)-fc,R 2+ C 2 R 2 3 + ^-lsinh(R a 4) 1 

rii R 2 [V J V H 2 J J 


(1.94) 
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Again the coefficients must be found from 


Q _ B 1 A ^ B g A 12 

^ A A A A 


^11 ^22 ^21 ^12 


and C 2 = B i A 2 i + B 2 A n 
”11 ^22 — ^21 ^12 


where now 


A n = cosh(R 1 L/2a)-cosh(R 2 L/2a) 

A 12 = R 1 2 cosh(R 1 L/2a)-R 2 2 cosh(R 2 L/2a) 

A 21 = R 1 sinh(R 1 L/2a)-R 2 sinh(R 2 L/2a) 

A 22 = R 1 3 sinh(R 1 L/2a)-R 2 3 sinh(R 2 L/2a) 

B i = cosh(R 1 L/2a)+^ 2 cosh(R 2 L/2a) 

B 2 = -^-sinh(R 1 L/2a)+^-sinh(R 2 L/2a) 

Rf R 2 


Having found the stream function as a function of E, and ri, the u and v may be 
found by 


u(5.tl) = - 


1 wm 

a 3ri 


and v(5,ti) = 


a a^ 


1.7.1 No Magnetic Field 

To test the validity of this method for obtaining an approximate solution to the 
biharmonic equation, the method will first be applied to two cases involving 
viscous dissipation terms; one in which the aspect ratio is unity, and the second 
in which the aspect ratio is large. Putting the stream function given by Eq. (1 .85) 
into Eq. (1 .24), we can identify b 2 = 8 and c 4 = 24. Therefore, Ri = 2. 1 09 +i 
0.6704 and R 2 = 2.109 -i 0.6704. 

Case 1 - Square cell with differentially heated sides. 

W=L=10cm, AT=10°C, |3=10' 5 /deg, g = 980 xl O' 6 , v = 3.46x1 0* 3 cm 2 /sec 
g II to L. 3p /3y = pAT/W, G = 0.00177. 
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The results are shown below: 


Approximate result 
*P 2.90 xIO -5 cm 2 /sec 

u 8.93 xl O' 6 cm/sec 

v 9.02 xl O' 6 cm/sec 


Numerical result 
3.58 xl O' 5 cm2/sec 
1.12 xl 0 5 cm/sec 
1.12 xl O' 5 cm/sec 


The approximate solution gives only fair results for the case of L = W since the 
assumption that the stream function has the same functional form over the entire 
length obviously breaks down in the corner regions. The maximum stream 
function is low by 19%, as is the maximum vertical velocity v. The maximum 
horizontal velocity u , is low by 20%. 

Case 2 - Vertical slot with differentially heated side walls. 


W = 2 cm, L = 10 cm, AT=10°C, p=10' 5 /deg, g = 980x10 ® 

v = 3.46x10- 3 cm 2 /sec, g II to L. 9p/9y = p AT/ W and G = 1.416x1 0 5 . The 
results are now: 


Approximate result 
'P 5.90 xIO -7 cm 2 /sec 

u 9.08 xl O' 7 cm/sec 

v 4.96 xl O' 7 cm/sec 


Numerical result 
5.99x1 0' 7 cm 2 /sec 
9. 1 6 xl 0 7 cm/sec 
5.60 xl 0 7 cm/sec 


The maximum value of the stream function is in excellent agreement as is the 
maximum vertical velocity, u. There is some disagreement in the maximum 
horizontal velocity, v, which can be attributed to the approximations used to 
obtain the 4>(x). However the approximation is only off by 1 1%. 

The accuracy of the approximation can be assessed by examining the residual 
field shown in Fig. 1 .22 . The normalized residuals were obtained by inserting 
the assumed solution back into Eq. (1.24), performing the indicated 
differentiations to evaluate the left-hand side at each point, and normalizing the 
result by dividing by the driving term. As can be seen, the only regions in which 
the approximation is in serious error is in the corners. This does not seem to 
significantly affect the overall accuracy based on comparison with numerical 
results (although the residuals obtained from numerical computations show a 
similar pattern), but it does explain why the approximate solution is more 
accurate for the larger aspect ratios. 


1.7.2 Effect of Strong Magnetic Fields 

The body force on an element of fluid in a magnetic field can be written 
F = cBx(Bxv) = o[(B*v)-(B»B)v]. 
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In two-dimensions this becomes 


F x = °[- B y 2v x+B x B y v y ] 

F v = a[-B x 2 v y+ B x B y v x ] 

Adding this contribution to the body force terms in the momentum equation, 

5p d 2 u 3 2 u D 2 

' a ; + ^ + ^ +p9, '' ( 2 +ovB ‘ B ’ =0 

~fy +,1 i>F + ^|^ + P 9 /- ovB > 2+c,uB « B » =° 

Introducing the stream function, cross differentiating, and subtracting to eliminate 
the pressure terms as was done in deriving Eq. (1.24), we get 


3x 2 dy 2 


3p 8p 


3x 2 x 


+ ^ B * 


2+2 S B > B « )=° 


Case 3 - Vertical slot with differentially heated side walls and a strong vertical 
magnetic field. Let the temperature distribution be given by ATx/W. Introducing 

scaled coordinates £ = x/a and q = y/a ( g taken along the negative x-axis) and 
using the Boussinesq approximation, the equation for the stream function 
becomes 

a 4 V a 4v F (3 AT g x a 4 2 3*V n 

a^ 4 afdrf an 4 w v d% (i.96) 

where the Hartmann number is defined as 

PoV 

Since the magnetic field does not interact with flow in the vertical direction, we 
can assume the core flow will have the same form as before, hence ^(^ri) can 

be written in the form <E>(£,)(ri 4 -2 r\ 2 -1). Putting this back into Eq. (1 .96) and 
performing the required differentiation, we identify b 2 = 8+Ha 2 and c 2 = 24 1/2 . If 
Ha 2 »1, both R 1 and R 2 will be real, but » 1 while R 2 « 1. This can cause 
some computational difficulty in determining the coefficients using the method 
prescribed above because of the large exponential values involved. Therefore 
some algebraic manipulation is required. Since R 2 « 1, Eq. (1.93) and (1.94) 
can be expanded about R 2 =0 to obtain 
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Fig. 1 .22 Plot of stream function residuals for Case 2. The approximate solution 
is reasonably good over the bulk of the field, but does break down in the corners. 
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+ 0(R 2 ) (1.97) 




_G 

R i 2 


(l* \ ( 4 \ 

+ C,~ C 1 + C 2 R 1 2 +-! f cosh(R^) 
V R i ) 


^- + 1 
2 


V * J 


®'(S) = 


$- 


Ci+C 2 R 1 2 + J 2 -jR 1 slnh(R 1 ^) 


0(R 2 ) 


(1.98) 


where G is given by 


Q _ P AT g x a 4 


Grv 


W v 16 


The coefficients Ct and C 2 are determined by requiring <E>(£ 0 ) =<i>'(£ 0 ) = 0 at = 
L /2 a. This yields 


and 


2R, 3 


*'(« = ^ 


R i 


R i(^o 2 - ) sinh(R 1 ) + 2 (cosh(R^) - cosh(R 1 ^ 0 )) 

sinh(R^ 0 ) 


sinh(R 1 £,)-^sinh(R 1 % 0 ) ' 
sinh(R 1 £ 0 ) 


(1.99) 


( 1 . 100 ) 


where R^ = Ha + V24/Ha + 0(Ha 2 ) =Ha, Ha » 1 . 

The horizontal velocity (in this case flow in the y-direction) is obtained from 

v(£,ti) = fx = iT° '( t 1 4 - 2 ^ 2 + 1 ) (1.101) 

and the vertical velocity (in this case flow in the x-direction) is 


u(^,r|) = = — <J>(t]-ti 3 ) 

3y a v ’ 


(1.102) 


For W = 2 cm, L 10 cm, G = 14.16, Ha =2860 1/2 , the results are : 


Approximate result 
0.0620 cm 2 /sec 
0.0953 cm/sec 
0.0248 cm/sec 


Numerical result 
0.0660 cm 2 /sec 
0.115 cm/sec 
0.0256 cm/sec 


'F 

0 

v 
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The stream function is shown in Fig. 1 .23. The effect of the vertical field is to 
round off the stream function along the centerline as can be seen better in the 

profile along r\ = 0 shown in Fig. 1 .24. Note that the stream function appears to 
approach the end walls at a finite slope, even though its derivative is zero at the 
wall. The effect of this is seen in the velocity profile shown in Fig. 1 .25 which 
shows the maximum velocity occurring almost adjacent to the end wall. 


From Eq. (1 .99), it may be seen that for » £ 0 , the stream function will have a 
maximum value of 


G | 

(L) 

i 2 _ Grv | 

rn 

| 2 _ Po 9 PAT W | 

f L T 

2 Ha 2 ' 

lw J 

1 32 Ha 2 1 

lw J 

! 8aB 2 1 

lw J 


The maximum vertical velocity is obtained from Eq. (1.102) 


(1.104) 


q _ 8 V = Gr v f L V 

3V3 a 6V3Ha 2 wiw; 


2p 0 gPAT f L y 
3V3aB 2 v Wj ' 


(1.105) 


Differentiating Eq. (1.100), we see that the maximum velocity horizontal velocity 
occurs at ti = 0 and when cosh(R 1 §) = sinh(R 1 ^ ) / R, ^ or ^ ) / R, £ 0 . 

The maximum horizontal velocity can be estimated by extrapolating the slope of 
the velocity curve to the end which gives 


« 1 G L Gr L v p 0 gBATf L ^ 

a 2 R z 2 2 8Ha 2 W W " 2aB 2 [wj‘ (1 ' 10£ 

This last result can also be obtained from the integral theorem. Taking a stream 
line just inside the walls, the driving term, £pg«ds = p 0 |3ATg x L . For large Ha, 
the viscous damping term is negligible compared to the magnetic damping term 
given by ^ aB 2 v • ds =aB 2 v 2 W . Equating and solving for v, gives Eq. 
(1.105). 

Note that both the horizontal and vertical velocity is independent of length scale 
for Ha » 1 . 


This formulation is also applicable to the case of a horizontal field applied to a 
horizontal slot with a linear thermal gradient along the length L. This 
configuration would roughly describe a Bridgman furnace with an axial field 
subjected to transverse accelerations. In this case 


q _ 9 y PATa 4 _ Grv W 
vL 16 L 
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Fig. 1 .23. Streamlines for Case 3 in which a strong magnetic field is applied 
vertically. The approximate result is on the left and the numerical result is on the 
right 
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Fig. 1 .24. Stream function along t| = 0 for Case 3. 
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and the above results would be multiplied by W/L which accounts for the reduced 
thermal gradient which is now AT/L instead of AT/W. In this case, u is the 
maximum horizontal velocity and v is the maximum vertical velocity. 

Case 4 - Narrow vertical slot with differentially heated side walls and a strong 
horizontal field magnetic field. The stream function equation, Eq. 1.95) simplifies 
to 


dx 4 


+2 


a 4 '? aV 


dx dy dy 4 J dy 




-^r-gx-a 


dy 


B„ = 0 


2 


Let the temperature distribution be given by ATx/W. Introducing scaled 

coordinates E, = x/a and r| = y/a, taking g along the negative x -axis and using the 
Boussinesq approximation as before, this becomes 


a 4 y , 2 a 4 v | a 4 y pat g x a 4 


a^ 4 afaq 2 at| 4 w v 


-Ha : 


a 2 ^ 

a 2 rj 


= 0 


(1.106) 


where the square of the Hartmann number is 


PoV 


For an infinite vertical slot, the ^-derivatives may be ignored and the r\ 
dependence of the stream function is given by 

|j-Ha 2 |j = G where G=-&*IS5l. 

an 4 ari 2 w v 


(1.107) 


With the usual B. C., \|/(±l) = \j/'(±l) = 0, the solution is 


V(q) = G 


cosh(Ha'n)-cosh(Ha) r| 2 -1 
Ha 3 sinh(Ha) 2Ha 2 


(1.108) 


Now an approximate solution for 'f'(£,q) is assumed to be of the form 
^(£.*1) = <K£)vOl) so that the equation for the stream function may be written 


an 4 J,. 0 a? an ! 


a> 2 av ^ 

„ =0 + aT *a? 


= G. 


(1.109) 


/ T ]=0 
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Evaluating the derivatives of y at r\ = 0 and rearranging we get 




l-cosh(Ha) 1 
[Ha 3 sinh(Ha) + 2 Ha 2 J 


+ 2 <j>" 


[Hasinh(Ha) Ha 2 


+[ 1]<|)=1 ( 1 . 110 ) 


If Ha »1, this equation can be simplified by retaining only the highest order 
terms in Ha, 


V 


2 Ha 2 


+ 2 <t>" 


Ha 2 


+[ 1]<|)=1 


or 


<|>""-4 <|) // + [2Ha 2 ](t)= 2Ha 2 . (1.111) 

This has the same form as Eq. (1.86) with b 2 = 4, c 4 = 2 Ha 2 , and G = 2 Ha 2 . 
Thus the solution can be written as the product of Eq. (1.108) and either Eq. 

(1.88) or Eq. (1.93) in which 

R _ V4 + 2(2Ha 2 ) V 4 ~2(2Ha 2 ) 

R _ V4 + 2(2Ha 2 ) V4-2(2Ha 2 ) 


For W = 2 cm, L 10 cm, G = 14.16, Ha =2860 1/2 , the comparison of the above 
approximate calculation and the numerical results are shown. 


Approximate result 
T 0.00248 cm 2 /sec 

u 0.00493 cm/sec 

v 0.00947 cm/sec 


Numerical result 
0.00249 cm 2 /sec 
0.00496 cm/sec 
0.00935 cm/sec 


Here the agreement is quite good, the difference being only a few percent. 


The stream function is shown in Fig. 1 .26. Note that the flow is rectalinear along 
most of the length of the ampoule and that the turning points are pushed very 
close to the ends. Since all of the change in the stream function takes place over 
a very short distance in the x-direction, it is possible to the much simpler semi- 
infinite formulation to obtain the same results. From Eq. (1.92), we can write the 
^-dependent part of the stream function as 


<K£)= 1- ( 7T s in(P^)+cos(pl;) 
vP 
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Fig. 1.26 Streamlines for Case 4 in which a strong horizontal field is applied to a 
vertical slot with differentially heated side walls. The approximate result is on the 
left and the numerical result is on the right 
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where a = Re(R1) and p = lm(R1). When £ > 1/a, this function quickly 
approaches 1 . The magnitude of the stream function then is determined by yfa) 
given by Eq. (1.108). Since this function has a maximum at r| = 0, the maximum 
value of is 


'f' = G 


l-cosh(Ha) 1 
Ha 3 sinh(Ha) + 2Ha 2 


Ha»1 


G 

*2 Ha 2 


Grv _ p 0 gpATW 
32 Ha 2 8aB 2 ' 


Actually, the peak value of 'F is slightly higher than this because of the small 
over-shoot when the pre-exponential factor goes negative before it is killed by the 
exponential term. This is reflected in the slight bulging of the stream function at 
the turning point. This effect is also seen in the numerical computations. 

The horizontal velocity along rj = 0 is given by 


v (5.o) = - 


ay 

dx 


G l-cosh(Ha) 1 
a Ha 3 sinh(Ha) + 2Ha 2 


P 2 sin(b^)e~ <,s . 


(1.113) 


For Ha » 1 , a * p = Ha /2 /2 1/4 and the last term reaches a maximum value of 
0.622 a at a £ = 0.7 . Therefore, the maximum v near the ends is 


v = 


G 

2aHa : 


-0.622a = 


Ha W 


0.523 =0.0322 


Grv 

Ha 575 ' 


( 

= 0.0925 p 0 p AT 

v 


PqW 2 

vo 3 B 6 


\ 1/4 

) 


(1.114) 


The maximum vertical velocity along the side walls away from the ends is 
obtained by differentiating y(ri) to obtain 


G sinh(Haii) 
aHa 2 sinh(Ha) ^ 


(1.115) 


Again the maximum will occur just inside ti = 1 and can be estimated by 
extrapolating the slope of Eq. (1 .1 15) to T] = 1 which yields 


u _ G _ Grv _ p 0 gPAT 
aHa 2 8Ha 2 2aB 2 


(1.116) 
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This result can also be easily obtained using the integral theorem by equating the 
driving term, p 0 g x (3 AT L to the magnetic damping term, aB* u(2L) and solving 
for u scales as W 1/2 . 


1.8 SUMMARY 

Methods for approximating buoyancy-driven flows in various geometries have 
been developed which are applicable over a wide range of Gr, Pr, and aspect 

ratio; limited only by the requirement that the Re be sufficiently low (typically Re < 

3000) so that the flows remain laminar. The most general result from Eq. (1.53) 
can be expressed for L < 2 W in the form 


16 y 

3V3 W’ 


~ J_ 

V " 3V3 8 e „ 


8 y ( 32 4* 2W 'j 

3V3W^27V3 k + L j 


where 


(1.117) 


4* 


/ « 

16 4 # 

K 

- 

[27V3 K 

4> 

r 

( 16 4* 

K 


[27V3 K 



\ 

+ 

J 


U J 


+ T_L_ = GrPr 
k W 384 

+ n = GrPr 
k W ~ 384 


for Pr > 1 


for Pr <1. 


In these expressions, L always refers to the horizontal dimension and W to the 
vertical dimension. The Gr is defined as g p AT W 3 /v 2 . 

For L > 2 W, these become 


u = 


16 4* 


3V3W’ V 3V35 e „ 3V3W 


84* 1 


8 4 / 


32 4* 

27 V3 k 


+ 1 


(1.118) 


where 


K 


r 16 y >3 

27 V3 k 


+ 'i 


4^J__ = GrPr 
+ k W “ 384 


for Pr > 1 


4 > 

K 


16 4 / 

27 V3 k 


r, ^ 


J 




-|3 


+ 4 / L _ GrPr 
k W ~ 384 


for Pr < 1 . 


Note that for large Gr Pr, these reduce to 
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H* = ^(6V3) 1M (GrPr) 1/4 K = 0.5049(6V3) 1/4 (GrPr) 1/4 k 


u = 


16 


2^3 


J/4 


GrPr 


w = 01943 


G^V /4 _v_ 
Pr 3 J W 


\V2 


v = ( 6 ^) (GrPr) 1/2 — = 0,2686 1 — 1 — . 
12 1 ’ W IPrJ W 


^ \ 1/2 

Gm v 


(1.119) 

( 1 . 120 ) 

( 1 . 121 ) 


This last equation has the form of the scaling law suggested by Ostrach for Gr 
>1, Pr >1, but should apply also for Pr <1 provided Gr Pr »1. The scaling laws 
proposed by Ostrach, Eq. (1 .2) and (1 .3), actually apply to flow in the vicinity of a 
vertical heated wall in a semi-infinite fluid which is different from the case of flow 
in a box with differentially heated end walls in which the streamlines are closed. 
An approximate solution to the vertical heated wall was obtained by Goldstein 
[12]and can be written as 


. 5.17 

v = 

3 


Gr 


J/2 


v 


pr+20/21; 


W 


( 1 . 122 ) 


which would scale as Eq. (1.2) in the limit of small Pr. That this scaling does not 
apply to the flow in a completely filled box is seen in Fig. 1.16 which plots log of 
Re 2 vs. log of Pr for constant Gr. Although the data for small Pr do not exactly 
fall 

on the 1/Pr slope or on the line predicted by Eq. (1.53), they certainly do not 
scale as 1/(Pr+20/21). 

For Gr Pr < 0(1 0 2 ), these reduce to 


Y = 


u = 


J_ 

384 

V 

16 'P 


W/L 
1 + (W/L) 4 

GrPr k 




v = 


16 W GrPr k 


and for L > 2 W they become 
GrPricf W/L 


y = 


384 


GrPrx 

/ 

W/2 < L < 2 W 

f 1 1 

Gr vj 

< 1 

U+w 4 /l 4 J 

~ 72V3 L 1 

J + W 4 /L\ 

( W/L > 

1 Gr v 

( W/L 

V 1 + W 4 / L 4 j 
\ 

‘ “ 72V3 L 

ll + W 4 /L 4 


1 + W/8L 


(1.123) 


(1.124) 


(1.125) 


(1.126) 
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- _ 16 T GrPr k 

f 1 1 

Gr v 

f 1 1 

U 3V3 W 72 V3 L 

[1 + W/8L J 

~ 72V3 L 

^1 + W/8L J 


8 T GrPr k 

f 1 ] 

Gr v 

r i ^ 

3-V3 W 144V3 L 1 

[l + W/8L J 

1 ~ 144V3 L 

^1 + W/8L J 


(1.127) 

(1.128) 


Note that these equations are slightly different from those derived earlier for the 
case of small Ra (Eq. (1.30) and (1.34)). This is because in the previous 
derivation the peak value of the stream function given by Eq. (1.25) was seen to 
give satisfactory results for all values of W/L and its shape was modified to 
account for the turning points being nearer the end for L > 2 W. The stream 
function used in the above set of equations was derived from the integral 
theorem by equating the viscous drag around the system to the driving force. 
The difference is small, however; amounting to a maximum of only 3% for L = 3 

W. In both cases, the 0 converges to the value obtained for maximum flow in a 
one-dimensional horizontal slot as L » W. 

Robertson and Spradley [6] numerically computed the maximum velocity 
generated by a horizontal temperature gradient in 2-dimensional chambers 
having a variety of shapes. They expressed their results in terms of a 
dimensionless velocity v* defined as 

* . 128v 

v = v - 

gPATd 2 

where d is the diameter of a circle of equivalent area. For a vertical slot, L = W/2 
and Ra = 1000, they obtain v* = 0.39. This would correspond to 

v = 0 30 9P AT 4W ' = 9PATW 2 
128 v n 2 515.54v ' 

Eq. (1.125) gives 

. _ Grv 2 2 gpATW 2 
72-V3 W17 " 530.01V ' 

For a horizontal slot, L = 2 W and Ra = 1000, they obtain v* = 0.18. This would 
correspond to 

v ■ 0 .i 9 9MLf 2 w* - 

128v Jl 264.55V ' 

Eq. (1.124) gives 
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0 _ Gr v 1 _ gpATW 2 

72V3 L (1 + 1/16) 265-OOv ' 

An upper limit for the time constant for transient flows was obtained by integrating 
around the streamline corresponding to the maximum steady state flow (Eq. 

(1 .39)). For the case of small Ra, this becomes 


w 2 

1 + (W/L) 2 

L 2 

1 + (L/ W) 2 " 

36 v 

1 + (W/L) 4 _ 

36 v 

1 + (L/W) 4 _ 


(1.129) 


Note that this function is symmetric with the exchange of L and W. The function 
reduces to W 2 /36 v for L = W and for L » W and has a maximum value of 1 .2 

W 2 /36 v for L = 1.6 W. The behavior is the same for L < W with the role of L and 
W exchanged. 

Robertson and Spradley also computed the response times for their various 
configurations which they expressed as a dimensionless time x* = t v/ R 2 , R 
being the radius of a circle of equivalent area. They obtain x* =.04 for the 

horizontal slot with L = 2 W, and x* =.03 for the vertical slot with W = 2 L. These 
would correspond to 


L 2 

X 52. 36 v ’ 


W = 2L 


and 


W 2 

X ” 39.27 v’ 


L = 2 W. 


These compare favorably with Eq. (1.129) which yields 


w2 L 2 
X “ 30.6v ~ 30.6 v’ 


L = 2 W and W = 2L. 


Eq. (1.129) assumed the maximum velocity was always at E, = 1/3 1/2 , as was 
discussed previously. Taking a more realistic value of E, to account for the fact 

that the maximum velocity occurs closer to x £ =1 would improve the 
agreement. The fact that the model fails to predict a difference between the time 
constants for the horizontal and vertical slot is probably due to the assumption 
made in deriving Eq. (1 .39) that the positions of maximum flow, L' and W\ stood 
in the same ratio as L and W. Even so, Eq. (1 .129) is accurate to within a factor 
of 2, which should be sufficiently accurate for order of magnitude estimates of the 
time constants. 
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The response to periodic accelerations may be approximated by rewriting Eq. 
(1.18) as 


v = v 


(1 + icox) 


= V- 


(1 + C0 2 t 2 ) 1/2 


(1.130) 


where |v| is the amplitude of the velocity fluctuations produced by a periodic 
acceleration whose magnitude would produce a steady state velocity v. The 

accuracy of this expression is also limited by the accuracy of x, but should be 
within a factor of 2. 


The effects of stabilizing and destabilizing axial gradients were considered. It 
was shown that the effect of such gradients on the flows induced by transverse 
accelerations could be represented empirically by 

u(Ra(x)) 500 

u 0 ~500 + Ra(x) C' 131 

where u(Ra(x)) is the maximum flow velocity with the stabilizing (or destabilizing) 
axial gradient, and u 0 is the maximum velocity without the axial gradient. The 
Raleigh number here is defined as 


Ra(x) = 


9x ft AT a 4 

KVX 


for dilute systems stabilized by thermal gradients, and 


Ra(x) = 


g x (Ap/p)a 4 

DvX 


for non-dilute systems stabilized by solutal gradients. In these definitions, X 
refers to the e-folding length of the axial gradient and Ra >0 is taken to be a 
stabilizing gradient. From this, it is seen that microgravity experiments involving 
low Pr fluids are relatively insensitive to the effects of axial gradients since IRal 
will generally be « 500. Non-dilute systems, especially those with large Sc, 
could benefit from a stabilizing gradient, and, consequentially, could be seriously 
affected by a destabilizing gradient. 

Finally, the effects of magnetic fields was assessed. It was shown that the 
application of a magnetic field can dramatically lower the convective flow velocity, 
even in a microgravity experiment where the flows are already extremely low. 
Axial fields are less effective than transverse fields for reducing flows, but have 
the advantage of being symmetric and also tend to keep the axial flows further 
from the solidification interface which lessens the radial segregation. The 
maximum horizontal flow velocity in a horizontal channel with a linear axial 
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temperature gradient in the presence of a strong axial magnetic field (Ha » 1) is 
given by 


u- Gr L v _ 2 Pq g ft AT f L > 
6V3Ha 2 WW 3V3oB 2 lwJ 


(1.132) 


Note that the flow velocity is independent of length scale. The ratio of velocities 
with and without the magnetic field can be written as 


Gr L v 
6V3Ha 2 W W 
u 0 Gr v 

72^3 L 

Having obtained the flow is only half the battle for most applications. There is still 
the question of the transport. As stated previously, the heat and mass transport 
can be characterized by the thermal and solutal Peclet numbers which are useful 
to estimate whether transport is dominated by convection or conduction, but 
there is no simple way to directly relate, for example, the solutal Peclet number 
to the degree of radial segregation in a solidifying ingot to the solutal Peclet 
number because of the complex interplay of a number of length scales including 
the diffusion length, the ampoule radius, and the turning point of the flow. This 
problem is the subject of the next two parts of this report. 


12 

Ha 2 


W 


Ha 2 »1 


(1.133) 
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SECTION 2 


RADIAL SEGREGATION IN BRIDGMAN GROWTH 

FROM 

LOW-LEVEL AXIAL ACCELERATIONS 
2.1 INTRODUCTION 

The importance of being able to control buoyancy-driven convective flows in 
order to grow compositionally homogeneous crystals from the melt has been 
recognized for some time and has motivated a number of experiments using the 
reduced environment of orbiting spacecraft to minimize such flows [1]. However, 
it has been found that even the reduction of the effective acceleration of gravity by 
some 4 to 6 orders of magnitude is not necessarily sufficient to eliminate 
compositional inhomogeneities in the growth of such crystals and the very slow 
flows driven by the residual accelerations inherent in a low Earth orbit spacecraft 
can be significant in crystal growth. 

A number of attempts have been made to estimate the magnitude and effects of 
such flows by scaling analysis. Camel and Favier identified various convective 
regimes based on the relative magnitudes of the growth Peclet number versus 
the product of the Grashoff and Schmidt numbers [2]. Langbein and Tiby 
developed criteria for allowable accelerations for a number of fluid and 
solidification processes based on scaling analyses [3]. Alexander and 
Rosenberger compare the results of the various scaling analyses to numerical 
simulations and find that such scaling arguments can at best only provide a rough 
order of magnitude estimate because of their neglect of the multi-dimensionality 
of the physical growth process [4]. 

Chang and Brown were the first to actually compute the buoyancy-driven 
convective flows in vertical Bridgman crystal growth of a dilute binary system and 
predict their effect on melt interface shape and dopant distribution [5]. Since then 
a number of researchers have utilized various computational fluid dynamics 
(CFD) techniques to extend this work to non-dilute systems, to non-axisym metric 
systems, and even to systems subjected to time-varying accelerations [6-9]. 

These studies have provided much useful information on the role convection 
plays in such growth processes, and have been particularly valuable in the 
planning of experiments to be conducted in the reduced gravity environment of 
an orbiting spacecraft in order to control or eliminate some of the unwanted 
effects of buoyancy-driven convection. Perhaps the most surprising result of 
these studies is the extent that one must reduce the effective acceleration level in 
order to actually achieve diffusion limited growth, even for systems of modest 
dimensions. 

One difficulty with the CFD approach is that a calculation only addresses a 
particular case and gives little insight into the role of the various parameters 
involved such as growth rate, the thermal profile, and the thermophysical 
properties of the candidate materials. An approximate analytical solution that is a 
function of these parameters would also be desirable to elucidate the complex 



interplay of the multiple length scales involved and to allow various sensitivity 
analyses to be performed without the expenditure of large amounts of computer 
time. 

In this section approximate analytical solutions are sought that capture the 
essential elements of the more definitive CFD models, but still remain 
mathematically tractable. In this first attempt, the analysis will be restricted to 
plane-front solidification of a dilute system in the vertical, thermally stable 
configuration (acceleration vector aligned along the furnace axis antiparallel to 
the thermal gradient) [10]. These approximate solutions may be compared with 
the CFD solutions at specific points in parameter space in order to determine the 
validity of some of the approximate methods that were required to obtain the 
solutions. They then may be used to extend the modelling over a much broader 
range of parameter space than could be reasonably carried out using CFD 
techniques. 


2.2 ANALYSIS 

The full set of coupled partial differential equations that govern the fluid flow, heat 
transport, and mass transport in even the simplest Bridgman crystal growth 
system indeed presents a formidable challenge. However, for many problems of 
practical interest it is possible to make a number of simplifying assumptions that 
allow closed form approximate solutions that are accurate to well within an order 
of magnitude that can provide insight to the roles of the various parameters. 

For example, since many of the materials of interest have low Prandtl numbers, 
heat is transported by conduction for low to modest convective flows in research 
scale growth systems. For dilute systems, the concentration field does not directly 
influence the flow field. This decoupling allows the flows to be estimated from an 
invariant thermal field, and the concentration field to be computed from the 
resulting flow field. Furthermore, since in Bridgman growth, the growth ampoule 
usually has a large aspect (length to diameter) ratio, a good estimate of the axial 
component of the convective flows can be made using a one-dimensional model 
which has a simple closed form solution. 

The solution of the mass transport equation is more difficult in that it is necessary 
to solve this partial differential equation in two dimensions in order to obtain the 
radial dopant distribution. For weak convective flows, however, it is possible to 
use perturbation theory to simplify the mass transport equation so that a solution 
can be obtained by separation of variables. Also, this analysis will be restricted to 
plane-front solidification to simplify the boundary conditions at the melt-solid 
interface. As was shown by Corriel and Sekerka [11], it is essential to arrange the 
thermal profile to achieve nearly a flat interface if radial segregation is to be 
minimized. 

In order to estimate the flows in a vertical Bridgman crystal growth system, 
consider the one-dimensional flow in a vertical cylindrical ampoule with a radial 
temperature distribution given by 
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T(r) = T 0 + (AT) (r/a) 2 


( 2 . 1 ) 

where To is the temperature at r = 0 and (AT) is the average radial temperature 
difference. 

The momentum equation is 

P + p(v- V)v = -Vp + pV 2 v + pg . 

at (2.2) 

In cylindrical coordinates the Laplacian operator for the z-component is 

y2y _ 1 aw ^ 3 2 w 

r 3r a r 2 ' 


Following the example given in [13], the inertial term, which is second order for 
the small velocities considered here, is neglected, the density is expanded about 
some reference temperature Ti, i.e., 


P = Pi + ( T -Ti) = pi -pi P [(T 0 - Ti ) + (AT)(r / af 


(2.3) 


where p the thermal expansion coefficient, and the pressure gradient is assumed 

to be due solely to the weight of the fluid, ie. Vp = - g pi where g is taken as 
positive downward. For steady acceleration, eq. (2.2) becomes 


3 2 w i 9w 
3r 2 +r 3r 


g P (at) 


(r/a) 2 


(Tq-Ti) 

(at) 


= 0 


(2.4) 


The boundary conditions are: 

w(a) = 0 (no slip at the walls). 


f 


w(r) r dr = 0 


and symmetry about r = 0. 
the quantity 


(no net flow). 

The conservation of flow boundary condition requires 
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and the steady state solution becomes 


w(r) = _ 9 . P ( AT ) a2 [-, _4( r / a )2 + 3 [r/af] . 

48 v (2.5) 

From this solution it may be seen that the maximum convective velocity will occur 
at r = 0 and may be written 

W = g a 2 

48 v (2.6) 

This could also be written as Re = Gr/48, from which it may be seen that simple 
scaling based on dimensional analysis over-predicts the maximum velocity by 
considerably more than an order of magnitude. 

Let us now consider the effect of the buoyancy-driven flows on the radial 
segregation in the sample. The mass transport is governed by, 


,,dC , 
u dF + (w 


/ ) dC n i d 
g ' dz D r dr 


dC 

dr) 



dZ 2 


(2.7) 


where v g is the growth velocity and u and w are the buoyancy-driven velocities in 
the reference system moving with the growth velocity. 

For plane-front solidification, the boundary conditions are given by 

c(r,°o) = Coo for all r 

and 

D (§) 0 + Vg(l " k )C( r f 0) = 0 f° r all r. 

In the absence of buoyancy-driven convective flows, the solution becomes 


C z = [(l - k) e' z/s + k] — where 5 = D/v g . 

k y 


( 2 . 8 ) 


Now consider the case of convective flows characterized by a velocity w « v g . A 
solutal Peclet number may be defined as P e = ^ 8 ID where 5 is the diffusion 



length of the solute boundary layer at the growth interface. Since 5 = D/v g , the 
Peclet number becomes just the ratio of the characteristic convective velocity to 

the growth velocity, or P e = w / v g . The velocities u and w can be expressed as 


u = v g P e u(x , r) 

and 

w = v g P e w(r,z) 


(2.9) 

( 2 . 10 ) 


where tf and W denote the dimensionless convective velocity components that 

have been scaled by the Peclet number. Also let the concentration field be given 
by 


0(r,z) = C z + P e CM). 

Putting these back into eq. (2.7) yields 

P* V 9 | + P.v 9 W 3 A + Pe * Vg * * . 


( 2 . 11 ) 


dz 


dz 


- P.a“ + P . D *? + D 


3 2 C 2 . „ „ d 2 C 


r 3r ” " 3r 2 dz 2 

and the boundary condition becomes 


+ Pe D 


dz 2 


( 2 . 12 ) 


lac, 

^ )z=Q 


9C 

ldz/z=0 


Dl^^l +p e D|^ + v g (l - kjC^o) + P e v g (l - k)c(r,0)= 0 for all r. 


The zeroth order terms in P e in eq. (2.12) are the terms in the unperturbed 
equation that is satisfied by C z . Dropping the terms of second order in P e , the first 

order perturbation in composition, C must satisfy 


^.v a w^ + D^ +D ^c +D a!c = 0 

dz 2 


9 dz Vg " 3z ' r 9r ' ~ 9 r 2 


(2.13) 


with the boundary condition 

d(||) q + v g (l -k)c(r,0)= 0 for all r. 
Dividing through by v g and identifying 6 = D / v g , 
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dC 

dz 


- W 


^ + &^ +5 ^ +8 ^ =0 

dz r 3r $ r 2 dz 2 


(2.14) 


Let us write 

w = w r w z 


(2.15) 


where W r and are respectively functions of r and z only. Let us also assume C 
can be written as 


C(r,z) = w r f(z) e‘ z/5 . 
The derivatives are 


§3^ = w r (f'(z>-M| e - 2 « 

dZ s 


(2.16) 


- 2 ^ r,z ^ = w r (f"(z) - 2 f ^ + e' z/6 . 

az 2 \ 8 6 2 / 


Putting these back into eq. (14) and dividing through by ffi r f(z) exp(-z/d), 

1_(f'( 2 ) - ML + J5-ff , ( z)- 2 f ( z ) + M| + 

(z)1 5 I f(z) 8 k f(z)\ 8 s 2 / 


1 

(«'(*) - 

f(z) 

+ 8 

1 1 dw 

Wr 

i r dr 


dr 2 


(2.17) 


We see that the variables separate provided 
1 J 1 dw r , d 2 W r | _ o2 

fiTU dr dr* / 

where s is the separation constant. 

Eq. (18), can be written 

r 2 d^w \ +r dw L+r 2 s 2^ r:=0 
dr 2 dr 

This is a form of Bessel's equation whose solution is 
w r = E J 0 ( r s ) + F Yq( r s ) . 


(2.18) 


(2.19) 


( 2 . 20 ) 
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Since Y 0 blows up at the origin, F must be set to zero. It is necessary to introduce 
a series solution to represent the velocity function, W r given by eq. (5), 


w 


©o 

r = I E nJo( r Sn) = [l - 4 (r / af + 3 (r/a) 4 ] 


n=1 


where the coefficients E n are computed using the orthogonality 
Bessel functions, i.e., 


E n = 


r f 

>n) Jo 


r (l - 4 (r/a ) 2 + 3 (r/a) 4 ) J 0 (r s n ) 


a 2 J^a s r 

which after performing the integration becomes 


E n = 


64 


1 


Ji(a s n ) [(a Sn ) 3 (a s n ) 5 . 


( 2 . 21 ) 
properties of the 

dr 


( 2 . 22 ) 


The separation constants s n are the zeros of Jo divided by the radius of the 
ampoule, ie. Si = 2.405/a, S2 = 5.520/a, S3 = 8.654/a, S4 =1 1 .792/a, etc. It was 

found that 4 terms were sufficient to give a reasonable representation of W r . The 
coefficients are Et = 0.332, e 2 = 0.898, e 3 = -0.335, and e 4 = 0.161. 


Now f(z) can be found by solving the ordinary differential equation, 


fn (z) - 


fnjz) 

5 


2 „ , (1 - k) C 

Sn 2 fc) =- L — ' — W Z 
s 2 k 


(2.23) 


In order to proceed, we must somehow specify $f z which, of course, could not be 
determined from the 1 -D solution of the velocity field. There are several 
conditions that this function must satisfy in order to represent the actual flow field. 

For example, <JJ Z should be zero at z = 0, peak at unity at some specified distance 
from the growth interface, and then eventually die out for large z. The behavior of 

<ff z at large z does not significantly influence the solution near the interface which 

is the region of interest. The distance at which <fJ z peaks is the point at which the 
axial flows begin to be turned by the pressure gradient resulting from the flow 
encountering the growth interface and will generally be near the point of 
maximum radial thermal gradient. Further, since u must be zero at z = 0 for all r, 
from the equation of continuity, 
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1 d(ru) dw 

r dr dz 


= 0 


it is seen that the derivative of Qf z as well W z must vanish at z = 0. Lacking a 
solution to the two-dimensional momentum equation a function with these 
properties will be arbitrarily assumed, ie. 



Here X represents the boundary layer thickness of the horizontal or radial flow 
across the growth interface [14]. The coefficient of eq. (24) is chosen to normalize 
the expression to unity at 2X where it achieves its maximum value. This will 

generally occur where the point at which the radial temperature gradient reaches 
its maximum value. 


For convenience, we write eq. (2.23) as 

f n (z) - a f n (z) - s n 2 f(z) = G z 2 e^ z 
where a = 1/8, b = VX, and 

G = - (l - k) a 2 p 2 ei — . 

4 k 


This ordinary differential equation has a solution given by 

f n(z) = An z 2 e-P z + B n z e‘P z +Cn e'P z + D n e fnZ (2.26) 

where 

A | _ G _ B _ 2 G (2 + g| 2 G (3 p 2 + 3 gp + g 2 + s n 2 ) 

p 2 + ap-s„2 " (p 2 1 a p . Sn2 f ' " (p 2 + a p-s„ 2 f 


and 


r n =|(l - Vi +4 s n 2 /a 2 ) . (227) 

The coefficients D n are arbitrary constants to be determined from the boundary 
conditions. (Since the perturbation must die out for large z, the other arbitrary set 
of arbitrary constants, i.e. the coefficients of the exponential of the positive values 
of r n , must be set to zero.) 
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The perturbed solution can now be written as 


~ 4 

qr , z) = X En Jo(r s n ) [<A n z 2 + B n z + C n ) + D n e«] e' 2/5 . 

n=1 


(2.28 


The arbitrary constants are found from the boundary condition at z = 0, ie., 

fac) 


v g (l - k)c(r,0) = 0 for all r. 


(2.29) 


Dividing through by D W r and identifying D/v g = 6 = 1/a, this becomes 


4 

I 

n=1 


X En Jo(r s n ) [Bp - (p + a)C n + (r n - a )D n ] + a(l-k) £ E n Jo(r s n )(C n + D n ) = 0 

n=1 


Collecting terms, 


~T 

^ E n Jo(r s n ) [b r - 13 C n + r n D n - k a (Cn + D n )] = 0 

n=1 


(2.30) 


In order for this to hold for all r, the bracket term must vanish for each n. This 
requires 


n _ *^n + (p + k a) C n 

u n 

r n - k a 


(2.31) 


We can now calculate C(r,z) from eq. (16) and in turn determine the 
concentration field using eq. (11). 


2.3 DISCUSSION 

Having obtained an analytic solution, even with a restricted range of applicability, 
it is now possible to make a number of general conclusions about the interplay of 
the various materials and processing parameters and their effect on radial 
segregation. 

To begin, since the maximum or characteristic convective velocity depends 
directly on the acceleration in the range over which the first order perturbation 
theory is valid, the perturbation to the concentration field, and hence the degree 
of radial segregation, is linearly related to the applied acceleration. 
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The magnitude of the perturbation at the growth interface depends on the value of 
the functions f n (z) evaluated at z = 0. From eq. (26) 


fn(0) — C n + D n 

which can be written using eq. (2.27) and (2.31) 


(2.32) 


fn(0) 


2 G 

P 3 + r n 

3 p 2 + 3 a P + a 2 + s n 2 ) + s n 2 (3 p + a) 

r pi - ka 

( 

P 2 +a p - s n 2 

3 


(2.33) 


This rather complicated expression involves the interplay of the three inverse 
length scales, ie. a, the inverse of the diffusion length; p, the inverse of the 
boundary layer of the horizontal flow across the growth interface, and s n , the nth 
zero of Jo divided by the radius of the sample. To gain some insight into the roles 
of these various length scales, it is instructive to consider several limiting cases. If 

a « s n , r n ~ -s n and eq. (2.33) becomes 


fn(0)-- 2 - 


_G 

Sn 


P 3 - S n 3 + S n 2 (3 p + a) - S n 1 

3 p 2 + 3 a p + a 2 

ii 

1 

p 2 + a p - s. 

2 

3 


(2.34) 


Recall from eq. (2.25) that G is given by 

G = -(l -k)a 2 p 2 — . 

4 k 


Using this, eq. (2.34) becomes in the limit of small a (small diameter ampoules in 
which the diffusion length 8 » a), 


fn(0)-^(l -k) 
2 k 


a 2 p 2 / |3 - s„ 


P - s, 


a « s n . 


n / 


(2.35) 


If p is large compared to s n , 


0»- — (1 -k) 
2 k 


or 


P S n 


or f n (o) ~ -3-^ , a « s n and p » s n 


If p is comparable to s n , this becomes 
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Coo e 2 

fn(0) = — — (l - k) — or f n (o) — ^ — , a » s n and p » a . 

2 k k p k 8 

From this analysis several conclusions can be drawn. For slow growth velocities 
in very small diameter ampoules, the diffusion length can be substantially greater 
than the ampoule radius, hence a < s n The n =1 term will then dominate 
because s n gets progressively larger with n. If the thermal profile is such that the 
radial gradients attain their maximum values at a distance comparable to or 
greater than a, the f n (0) term will go as the fourth power of the diameter and will 
be inversely proportional to both the square of the diffusion length and the square 
of the boundary layer thickness of the radial flow. Physically, this has the effect of 
lessening the radial concentration gradient by moving the radial flow farther from 
the solidification interface. Also in the limit of small diameters, the radial back- 
diffusion terms in eq. (2.14) become more effective in reducing the compositional 

perturbation driven by vb. Under these circumstances, the perturbation can be 
reduced by lowering the growth velocity. Even though this increases the Peclet 
number which multiplies f n (0), the function f n (0) goes as the inverse square of d 
which in turn is inversely proportional to the growth velocity. The net result is that 
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the radial segregation is directly proportional to growth velocity in this regime. 

The perturbation may also be reduced by moving the point at which the vertical 

flow becomes maximum farther from the growth interface, i.e. increase X. This can 
be accomplished by lengthening the adiabatic zone between the hot and cold 
zones of the Bridgman Stockbarger furnace. This also tends to flatten the 
isotherms in the vicinity of the growth interface and reduces the overall radial 
thermal gradient which drives the convective flows. 

For faster growth systems in larger diameter ampoules, the diffusion length may 
be small compared to the radius, hence a » s n . Thermal considerations 
generally require that the point of maximum radial gradient is comparable to or 
greater than the radius, thus X » 8 or a » p. Now as seen in eq. (2.37), the f n (0) 
no longer depends on the ampoule radius and each term will have the same 
magnitude. (The series is truncated in this case by the E n coefficients.) Also the 
magnitude of the perturbation is increased by the additional factor k in the 
denominator (assuming k<1 ). The f n (0) term will then be proportional to the 

square of the diffusion length 8 and inversely proportional to the square of the 
horizontal flow boundary layer thickness X. Physically, this may be understood 

from the consideration that if 8 « X, only a small portion of the concentration 
build-up at the solidification interface extends into the radial flow field. The 
perturbation may be lessened by further reducing the diffusion length by 
increasing the growth velocity, or by moving moving the radial flow field farther 
from the solidification interface by increasing the length of the adiabatic zone. 
Since the perturbation term is multiplied by the Peclet number, which is inversely 
proportional to the growth velocity, the radial segregation will now be proportional 
to the inverse cube of the growth velocity. It should also be remembered that 
even though the perturbation term is independent of a, the characteristic 
convective velocity, hence the Peclet number, and therefore the radial 
segregation, is proportional to the square of the radius. 


2.4 COMPARISON WITH CHANG AND BROWN'S RESULTS 

The very simple 1-D calculation for the velocity profile and the first-order 
perturbation calculation for the resulting concentrations may be checked against 
the detailed CFD computations carried out by Chang and Brown (C&B) for the 
flows associated with a dilute alloy system (Ga-doped Ge) in a vertical Bridgman 
configuration. Chang and Brown characterize their computations in terms of a 
Rayleigh number defined by 

Ra = 9 lie) L 3 (2.38) 

aj v 


where the length of the ampoule L = 2 cm, the radius a = 0.5 cm, the thermal 
diffusivity aj = 0.078 cm 2 /sec, the kinematic viscosity v = 0.0013 cm 2 /sec, and the 
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thermal expansivity p = 0.00025 /K. Values for g and Th - Tc were not specified, 
but can be related to the Rayleigh number by inserting these values into eq. 
(2.38), 


9 (Th - Tc) = 0.0507 Ra (cm K/sec 2 ). (2.39) 

In the growth of systems such as Ga-doped Ge a typical temperature difference of 
200 K is needed between the hot and cold zone to provide sufficient thermal 
gradient to stabilize the growth interface. Hence a Rayleigh number of 10 would 
correspond to an acceleration of 2.5 micro-gs. 

The quantity needed to evaluate eq. (2.6) is not the maximum temperature 
difference, but some form of average radial difference - (AT). It may be seen from 
the isotherms computed by C&B that the wall temperatures rise from Tc to Th in a 
nearly uniform manner over the length of the adiabatic zone, L a ; whereas the 
temperature along r = 0 rises from Tc to Th over roughly twice the length of the 
adiabatic zone. With these observations, a very rough estimate can be made of 
the average radial temperature gradients by 


J r 2 u 

[T(a,z)-T(0,z)] dz 

0 


Th-Tc 

8 


(2.40) 


C&B present their flow calculations in terms of stream function. The stream 
function y is defined such that the velocity vector can be obtained from 

w = - 7 - 3 ^- and u = 1 ^ . (2.41) 

r dr r dz 

The stream functions obtained by C&B were fit by even order polynomials to 
assure symmetry about r = 0 and the velocities were determined by 
differentiation, i.e., 

w(r) = 2 a 2 + 4 a 4 r 2 + 6 a 6 r 4 ....+ j aj ri ' 2 (2.42) 

where & 2 , a 4 , ... aj are the coefficients of the j ,h order terms in the polynomial fit 
for \| /. Chang and Brown state they used a value of 16 p.m/sec for their growth 
velocity v g ; however, this value seems inconsistent with both the stream function 
plots and the concentration fields they obtained. For the case of no convection, 
the growth or translational velocity will contribute v ? r 2 /2 to the stream function. 
The maximum dimensionless stream function obtained for this case (taken from 
their Fig. 2.4) was 0.0003. They scaled length by L and velocity by aj/L; 
therefore, their maximum dimensional stream function for this case would be 
0.000047 cm 3 /sec. At r = 0.5 cm, this would correspond to a translational 

velocity of 4 jiim/sec. 
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Chang and Brown chose a value of 0.00013 cm 2 /sec for their diffusion 
coefficient. For their stated value of 16 pm/sec for v g> the diffusion length d 
would be 0.08 cm. However, their concentration profiles indicate a diffusion 
length of approximately 0.2 cm which would correspond to a growth velocity of 
6.5 pm/sec. 

Because of these inconsistencies, it is difficult to make a direct comparison of 
flow velocities computed from eq. (2.5) with those corresponding to the stream 
functions presented by C&B. One approach would be to eliminate the scale 
factor by computing the ratio of convective velocities to translational velocity v g 
and then assume various translational velocities in making the comparisons. 
Table 2.1 shows the comparisons for the maximum convective flow velocities 
(downward flow at r = 0) after subtracting the translational velocity for the three 
values of v g in question. 


Table 2.1 


Ra 

Eq.(5) 

C&B 

C&B 

C&B 



vg=4 pm/sec 

Vg=6.5 pm/sec 

Vg=1 6 pm/sec 


w(0) (cm/sec) 

w(0) (cm/sec) 

w(0) (cm/sec) 

w(0) (cm/sec) 

100 

-0.00063 

-0.0005 

-0.0009 

-0.002 

1000 

-0.0063 

-0.005 

-0.009 

-0.02 

10000 

-0.063 

-0.05 

-0.09 

-0.2 

1000000 

-6.3 

-1.1 

-1.79 

-4.4 


The velocities estimated from eq. (2.5) fall between those computed from C&B's 
stream function for v g = 4 pm/sec and 6.5 pm/sec for Rayleigh numbers up to 
about 10 4 . For materials with low Prandtl numbers the heat is transported 
primarily by conduction at low flow velocities. As the Rayleigh number 
increases, the isotherms eventually become distorted by the convective flows 
and at Rayleigh numbers between 10 4 and 10 5 the heat begins to be 
transported more by convection. At around this point, the velocity will begin to 
increase as the square root of Ra, as can be seen by the departure of C&B's 
results from a linear relation with Rayleigh number at Ra = 10 6 , and the simple 
analytical approach used here is no longer valid . 

In order to compare the estimates of the radial segregation, the values for the 
segregation coefficient k, the diffusion coefficient D, and the growth velocity v g 
were chosen to match those used in Chang and Brown's calculation. C&B took 

k = 0.1 , D = 0.0001 3 cm 2 /sec, and v g = 1 6 pm/sec for However, as stated 
previously, this value for v g leads to an inconsistency in the diffusion length in 
the concentration profiles presented by C&B. Therefore calculations to be used 
in the comparison were performed with values of v g = 6.5 pm/sec as well as 16 
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pm/sec. 


The only other parameter needed for the first-order perturbation calculation is 
the distance 2A,from the interface at which the axial convective velocity attains 
its maximum value This will occur at the point where the radial temperature 
gradients are the largest, which usually occurs just above the adiabatic zone. In 
C&B's work the length of the adiabatic zone was taken to be 0.5 cm and the 
interface was approximately in the center of the adiabatic zone, hence the 

adiabatic zone length extended 0.25 cm into the melt. The length 2X was taken 
to be 20% above this point or 0.3 cm, thus X = 0.15 cm. 

With these inputs, the calculated parameters for Ra = 100 and k = 0.1 are given 
in Table II along with the radial segregation £ defined as the difference between 
the maximum and minimum concentration at the interface divided by the 
unperturbed interfacial concentration. 



Table 2.2 
v q . =6.5 pm/sec 

v fl =16 pm/sec 

vc (pm/sec) 

-6.348 

-6.348 

Pe 

-0.977 

-0.397 

8 (cm) 

0.200 

0.081 

f,(0) 

0.3479 

0.8828 

1.(0) 

0.0461 

0.1579 

1.(0) 

-0.0042 

-0.0184 

M0) 

0.00077 

0.0036 

C(0,0)/C«o 

0.1536 

0.4365 

C(0,0)/Coo 

8.50 

8.26 

e 

15.0% 

17.3% 


The perturbed concentration profile at the interface is shown in Fig. 2.1 and 
compared with the result obtained by C&B. Notice that the concentration 
profiles are quite similar except for the shift along the z-axis. Inherent in the first 
order perturbation approach is the assumption that the perturbation in 

concentration C could be written 
G(r,z) = w r f(z) e‘ z/5 . 

Since w r must be zero at r = a (no slip condition), the perturbation must vanish 
at the wall. 

The terms containing the radial velocity which transport the melt from the center 
of the ampoule to near the wall in the vicinity of the growth interface are second 
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Concentration Profiles at the Interface 



Figure 1. Comparison of the concentration profiles at the solidification interface 
computed from the analytical mode with those computed numerically by Chang 
and Brown. The shapes are quite similar; however, the results from the 
analytical model are displaced. This is the result of the neglect of the second 
order radial convective transport term in the first order perturbation calculation. 
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order and were therefore neglected in the first order calculation; thus allowing 
the variables to be separated. As may be seen in Fig. 2.2, these terms would 
have had the effect of removing some of the solute-poor melt from the center 
and transporting it radially to dilute the solute-rich melt near the wall, thus 
shifting the isoconcentration curves upward along the z-axis. 

The error caused by neglecting these terms is not serious for small values of Pe, 
but for the above case in which the Pe is nearly unity, the first-order perturbation 
model has clearly exceeded its range of applicability. However, the predicted 
values of radial segregation from the first order calculation, (15.0%, and 17.3% 
for the two growth velocities considered) are remarkably close to the value of 
14.6% estimated from Fig. 9 in Chang and Brown’s work. The reason the first- 
order calculation continues to give reasonable results is that the radial 
segregation is computed from the difference between the highest and lowest 
concentration at the interface, hence would not be significantly affected by a 
shift of the concentration profiles along the z-axis. 

Comparisons were also made with Chang and Brown's calculations of the 
effects of varying some of the material and processing parameters. In the 
following examples, the reference case is Rayleigh number Ra = 100, 
segregation coefficient k = 0.1, Schmidt number Sc = 10, sample radius a = 0.5 
cm, and adiabatic zone length L a . = 0.5 cm. The effects of varying the 
segregation coefficient k are shown in Table III. 


k 

Cc&B 

fable 2.3 

C(6.5um/sec) 

C(l6um/sec) 

0.1 

17% 

15.0% 

17.3% 

0.3 

10% 

9.4% 

8.2% 

0.5 

7% 

5.6% 

4.3% 

0.7 

3% 

2.9% 

2.0% 

0.9 

2% 

0.9% 

0.6% 


The Schmidt number was varied by varying the diffusion coefficient. The results 
are shown in Table 12.4. 


Sc 

Cc&B 

Table 2.4 

C(6.5M.m/sec) 

C( 1 6jim/sec) 

30 

36% 

47.4% 

18.4% 

20 

32% 

36.2% 

19.9% 

1 0 

16% 

15.0% 

17.3% 

5 

4% 

4.6% 

8.3% 


The effect of varying the length of the adiabatic zone was also investigated and 
the results are shown in Table 2.5. 
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Concentration Profiles 



Radial Distance (cm) 


Figure 2. Concentration profiles throughout the melt in the case considered by 
Chang and Brown. 


La(cm) 

2X(cm) 

Table 2.5 
£c&b 

C(6.5|im/sec) 

6jim/sec) 

0.5 

0.6 

20% 

15.0% 

17.3% 

1.0 

1.2 

6% 

8.7% 

7.7% 

1.5 

1.8 

2% 

4.6% 

3.2% 

2.0 

2.4 

<1% 

3.0% 

2.3% 


The stated Chang and Brown results were scaled from the figures in their paper 
and are probably not accurate to better than a few percent. Also, lengthening 
the adiabatic zone alters the thermal profile, reducing the average radial 
thermal gradient which drives the flows, resulting in smaller convective 
velocities. This was not taken into account in the above calculations which 
could account for the increasing discrepancies between this model and the 
results of the Chang and Brown model as the adiabatic zone is increased. 


2.5 COMPARISON WITH ALEXANDER, OUAZZANI, AND 
ROSENBERGER'S RESULTS 

A more extensive numerical study of the effects of low level residual 
acceleration on various microgravity experiments was recently completed by 
Alexander, Ouazzani, and Rosenberger (A.O.R.). Since they were more 
interested in non-axisymmetric cases as well as in time dependent cases, they 
solved the flow and mass transport equations in a two-dimensional rectangular 
geometry in order to simplify the problem and to keep the computer time within 
reasonable bounds. Several static three-dimensional cases were solved 
numerically as checks on the use of the two-dimensional model and the 
differences for the most cases were not found to be significant. They chose the 
same growth system as Chang and Brown and calculated the radial 
segregation for a variety of applied accelerations as well as for different sizes, 
temperature fields, and growth velocities. 

A comparison of the perturbation model with their results for cases in which the 
acceleration is perpendicular to the growth interface is shown in Table 2.6. 
Since A.O.R. assumed the the same thermal field as Chang and Brown, the 
same expression relating average radial temperature gradient to the hot zone 
and melting point temperatures, eq. (2.40), was also used in computing these 
comparisons. The other parameters were identical to those used in the Chang 
and Brown model. 

The agreement between the first order perturbation model and the results 
obtained by Alexander, Ouazzani, and Rosenberger is extremely close (within a 
few percent) for small Peclet numbers where the assumptions made in the first- 
order perturbation theory are strictly valid. For Peclet numbers approaching 
and even exceeding unity the differences are within a factor of two.. It is 
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interesting to note that the disagreement between the first-order analytical 
model and the computational results becomes more severe if the Peclet number 
is increased because of increasing g-level, which increases the actual 
convective velocity, rather than by decreasing growth velocity (e.g., compare the 

1 0‘ 4 g case for v g = 6.5 pm/sec with the 1 0' 5 g case with v g = 0.65 pm/sec). This 
is a result of the fact that the error introduced by dropping the radial convective 
transport term in the first-order model becomes more serious as the convective 
velocity increases. At the very low convective velocities, this transport would be 
reduced by the radial back-diffusion term. 


TABLE 2.6 

A.O.R.'s results This work 


g/go 

T H (K) 

Vq(ji m/sec) 

a(cm) 

C 

Pe 

c 

10-5 

1331 

6.5 

0.5 

6.4% 

-0.377 

5.8% 

io - 5 

1331 

0.65 

0.5 

0.95% 

-3.77 

0.80% 

5x10-6 

1331 

6.5 

0.5 

3.2% 

-0.188 

2.9% 

10- 4 

1346 

6.5 

0.5 

36% 

-4.34 

66% 

io - 5 

1346 

6.5 

0.5 

7.5% 

-0.434 

6.6% 

1.4x10-6 

1346 

6.5 

0.5 

2.0% 

-0.0607 

0.93% 

10-6 

1346 

6.5 

0.5 

0.7% 

-0.0434 

0.66% 

10-5 

1346 

3.25 

0.5 

4.6% 

-0.868 

4.1% 

10-5 

1346 

0.65 

0.5 

0.7% 

-4.34 

0.93% 

10-6 

1346 

3.25 

0.5 

0.4% 

-0.087 

0.41% 

10-6 

1346 

0.65 

0.5 

0.0% 

-.434 

0.09% 

IO’ 6 

1346 

6.5 

1.0 

3.8% 

-.174 

6.4% 


It is particularly gratifying to observe that the perturbation model responds to 
variations in the growth velocity in the same manner as the numerical model. 

On the one hand decreasing the growth velocity increases the Peclet number 
which increases the perturbation to the composition field. However, this also 
reduces the gradient in the diffusion field in front of the interface which lessens 
the compositional variations, thus offsetting the increase due to the larger Peclet 
number. The perturbation model appears to have correctly captured this 
feature. 

For the last case, in which the size was doubled, it was assumed that the entire 
system was scaled in proportion. This has the result of doubling the length of 

the adiabatic zone, hence A. However it is not clear whether this scaling 
actually preserved the thermal profile used in the the A. O. R. computations. 

This might account for the apparent difference between the two models which is 
larger than might be expected given the small value of the Peclet number. 
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2.6 RESPONSE TO TRANSIENT AND PERIODIC ACCELERATIONS 

The time-dependent equation for the z-component of velocity, eq. (2.4) may be 
written 


aw _ 3 2 w | i aw | g (3 (at) r , « _ j_ 

at ^ r 2 r a r V [ l 3 . 


(2.43) 


where the inertial term has been neglected, the pressure gradient was assumed 
to equal the rg term, and the (T 0 - T 1 )/(AT) term has been set to 1/3 to adjust the 
buoyancy term to provide no net flow. 

For a transient acceleration applied to a system initially at rest, ie. g is “switched 
on” at t = 0, the solution may be written 


w(r,t) = w(r) + w E n J 0 (s n r) e* s n 2 v * 

n=1 


(2.44) 


where w(r) is given by eq. 2.5. 

Notice that the time constant is given by 


tn “ . (2.45) 

s n 2 v v ' 

The most persistent term is n = 1 . The time constant associated with this term is 


Xn a? _ a 2 

2.4048 2 v 5.78 v 

which for a = 0.5 cm and v = 0.0013 cm 2 /sec, x = 33 sec. However, the initial 
response is much faster because of the higher order terms. As may be seen in 
Figure 3, the time required for the fluid to reach 1/e of its final velocity calculated 
from eq. (2.44) for the above case is only 12 sec. For comparison, the time 
constant estimated by Langbein and Tiby [4] was a 2 /n which would give 192 
sec. 


Griffin and Motakef (G&M) analyzed numerically the response of the melt in a 
Bridgman growth system to both transient and periodic accelerations [15]. For 
small Rayleigh numbers, they find the time required for the upper convective roll 

to reach 99% of its steady state value was 0.17 a 2 /v, which for the case at hand 
would correspond to 33.7 seconds (the 1/e time was not specified). Eq. (2.44) 
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would predict 117 sec to reach 99% of the steady state velocity. Since the flows 
in the upper convective roll closely correspond to the flows predicted from the 
one-dimensional model (G&M state that these flows are exclusively controlled 
by melt radius and kinematic viscosity and, to the first order, are independent of 
the longitudinal extent of the the convective cells), it is difficult to understand 
why the time response time could be significantly shorter than that predicted by 
the one-dimensional analytical model. 


For periodic accelerations, assume g(t) = $ e' 03 * where ^ is the amplitude of the 
acceleration. After the initial transient dies out, the velocity will be given by 

w(r,t) = Aw(r) e' 03 * where Aw(r) is the the complex amplitude of the velocity. Now 
eq. (2.43) becomes 

d 2 Aw + 1_ d Aw . i co Aw + 9 P (at) 
dr 2 r dr v v [' a ' 3 

It is convenient to express the solution as a Fourier-Bessel series 


(2.46) 


Aw(r) = w X En Jo(s n r) 

n=1 


(2.47) 


where w is the maximum steady state velocity given by eq (2.6) that would result 

from an acceleration amplitude c^, and E n ’ are the new complex Fourier-Bessel 
coefficients. This choice of representation of the solution automatically satisfies 
the “no-slip” boundary conditions at the walls. Taking derivatives, putting them 
back into eq. (2.46), and multiplying by a 2 , we get 


OO 


+ i Q 


En Jo(Sn r ) — 48 


n=1 


(sr-i 


(2.48) 


where the dimensionless frequency Q is w a 2 /v,and £ n are the zeros of Jo. 

Using the orthogonality properties of Jo, we can multiply by Jo(s m r) and extract 
the Fourier-Bessel coefficients; 


-84 

' 1 


2 \ 

1-iG/U 




\1+Q /^ m j 


(2.49) 


Note that Q/£ m 2 = arc m where x m was the time constant for the term in the 
transient solution. It may also be seen that for Q = 0 the solution becomes the 
steady state solution given by eq. (2.22) as E m ’ reduces to E m . 
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The amplitude of the velocity is 


Aw(r) = w {l [ Re(E m ) J 0 (s m r)] 2 + X [ I m(E m ) J 0 (s m r)] 2 } 


(2.50) 


The velocity amplitude at r = 0 is given by 


Aw(o) = w 


I 

Em 1 

2 

+X 

E m £2 

2 

i+n 

_U 2 (l +£2 2 /U 4 ). 





(2.51) 


If Q » t, m 2 , the imaginary contributions will dominate (eventually £ m 2 > Q, but by 
them the E m terms will no longer be significant) and the velocity amplitude at r = 
0 becomes 


Aw(o) w Y 

m=1 ^ 


>> 


(2.52) 


If we set Q = 0 in eq (2.48) for r = 0, we see that 
£ Em U 2 = - 16 

Therefore, for the higher frequency disturbances, the velocity amplitude at r = 0 
becomes using eq. (6) 

Aw(o) = w 1-2- = w 16_v_ _ Q » ( 2 . 53 ' 

^ co a 2 3 co 

Figure 2.4 shows how the velocity amplitude starts to diminish at around 0.01 
rad/sec and becomes inversely proportional to w at 1 0 rad/sec. For 

comparison, G&M find that for co » a 2 /v, 

Aw = S_v 5.9 w v 

0.17 co a 2 co a 2 


This has the same functional form as eq. 2.53, but differs by a factor of 2.7. One 
possible reason for this discrepancy is that G&M compared the amplitude of the 
fluctuating stream function with the value for steady accelerations, whereas this 
work compares the velocity amplitude at r = 0 with its value for steady 
accelerations. Even though the buoyancy driven convective velocity is 
maximum at r = 0 for steady accelerations, this is not necessarily the case for 
higher frequencies. Even so, the one-dimensional model seems to give close 
enough agreement to be useful for the purpose intended. 


94 



W(o,t) 

W(o,°o) 



Time (sec) 


Fig. 2.3. Build up of convective flow at r = 0 from a steady acceleration that is 
suddenly “turned-on” at t = 0. Even though the time constant for the dominant 
(most persistent) term is 33 sec, only 12 sec is required to reach 1/e of the final 
steady state value. 
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2.7 EFFECT OF PERIODIC ACCELERATIONS ON MASS 
TRANSFER 


The time dependent mass transfer equation is 


dC 

at 


+ (U ' Vg) 


ac 

ar 


+ w 


ac 

dz 


= D 


i ac a 2 c 

-I h 

r ar a r 2 


a 2 c 

az 2 


(2.54) 


Using the perturbation approach as before, the u and v velocity components are 
expanded in terms of the Peclet number as before, except now the Peclet 
number is assumed to be 


P G = AP e e icot = 


Aw(o) 


iicot 


(2.55) 


and the concentration is assumed to be expressed by 

C(r,z,t) = Cz(z) + AP e C(r,z) e icot . (2.£ 

Putting this back into the partial differential equation and dropping all but first 
order terms in APe as before, we obtain 


- i co C + v 


ac 

9 az 


-v g w 


aCz + Da6 D a^c 

az r ar dr 2 



dz 2 


= 0 


(2.57) 


The variables separate as before but eq. (2.23), the second order ordinary 
differential equation that must be solved for f n (z), now has an extra term, vis 


Mz) - 




+ ^|fn(z) =- 


(l -k) C» 
5 2 k 


(2.58) 


The solution for f n (z) will still given by eq. (26) if the s n 2 terms in A n , B n , C n , D n 
and r n are replaced by s n 2 +ico/D. 


The perturbed concentration field is then given by 


AP e C(r,z) =5S E n J o(Sn r)f (z) 

v g n 


(2.59) 


where E n ’ is given by eq. (2.49) above. 
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The fall-off of the perturbed concentration field at the interface with frequency is 
also shown in Figure 2.4. The controlling factor is the relative magnitude of s n 2 

and w/D. If co/D £ s n 2 (or co a 2 /D £ £ 1 ), then the concentration field cannot cannot 
change as rapidly as the applied acceleration and the amplitude of the 
oscillation will be diminished. As may seen this begins to occur at around 1CH 
Hz which corresponds to a value of co a 2 /D = 1 .28. The fall-off becomes more 
rapid with increasing frequency, and when the velocity amplitude also begins to 
fall-off with increasing frequency, the concentration perturbation fall as the 
inverse cube of frequency. 

This behavior may be understood by examining the various coefficients that 
make up the fn terms by replacing each s n 2 term with co/D. This yields 


- 2 G D 2 (2 (3 + a) 2 G D 2 r “ {l + i) 

co 2 CO 2 2 fi5 

- 12 GD 6,2 (3P+(1 +k)aKl - i) 

u n — » — 

0) 5/2 

»„(0) -* 2 G D 2 + i V2 G D 5/2 (3 p + (l + k)a) 

n rr.2 r.i 5 / 2 


In the limit of large co, the amplitude of the perturbed concentration field at r = 0 
on the interface becomes 


APe 6(0,0) = X E n f n (0) • 2 i G D 2 W £ En 

v 9 co 2 ^Vg 


(2.60) 


Using the definition of £2 and the previously determined value for the sum, this 
may be written 

APeC(0,0) -> 32iQ D 2 V w 

co 3 a 2 v g < 2 ' 6 

where G is given by eq. 2.25. 
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Log Frequency (Hz) 


Fig. 4. Response of the velocity and of the concentration at the interface to 
periodic accelerations. For very low frequencies (< 10* 4 Hz) the flow and 
concentration fields retain their steady state values. The velocity field will start 
to diminish at around 0.01 Hz and falls off as the inverse of frequency beyond 
0.1 Hz. The perturbation to the concentration field starts to diminish at around 
10' 4 Hz and falls off as the inverse cube of frequency above 0.1 Hz 
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The composition at the interface will fluctuate with a frequency w and the 
relative variation at r = 0, which in this case would be the radial segregation, is 
given by 


_ AP e C(0,0)l _ 8 e 2 (l -k) D 2 v w 


Cz(0) 


8 2 A 2 co 3 a 2 v n 


(2.62) 


Using the definition for d and eq. (2.6) for the maximum convective velocity, this 
may be written 


C - 


(l - k) p (at) g v 
6 X 2 to 3 


— , co » E,i 2 D/a 2 and to » ^i 2 v/a 2 


.(2.63) 


Note that this is independent of material parameters (except for expansion 
coefficient) and the only size parameter that appears is the length of the 
adiabatic zone. This is somewhat different from the order of magnitude analysis 
of Langbein and Tiby who estimated 

3C _ 9 A P/P w hich corresponds to C = ^ 

C L to 2 L to 2 

where L is some length scale. 


2.8 CONCLUSIONS 

It has been shown that a simple one-dimensional flow model can be used to 
obtain a reasonably accurate estimate of the maximum buoyancy-driven 
convective velocity of a dilute binary melt in a vertical Bridgman crystal growth 
system. Using this estimate, a first-order perturbation solution to the mass 
transport equation was obtained for convective flows that are small compared to 
the growth velocity. The radial segregation computed by this combined 
analytical model compares quite favorably with the numerical models 
developed by Chang and Brown and by Alexander, Ouazzani, and 
Rosenberger indicating that the essential physics has been captured by this 
model. This means that the model can be used with some confidence to predict 
the effects of varying material and/or processing parameters for experiments 
subject to small convective flows such as those produced by the residual 
acceleration on an orbiting spacecraft. It has also been found that even though 
the first order perturbation calculation of the concentration profiles is no longer 
valid as the convective velocities approach the growth velocity, the model 
continues to give reasonable values for radial segregation since the primary 
effect of the higher order terms is to shift these profiles along the vertical axis, 
thus the difference between their maximum and minimum value is not greatly 
affected until the convective flows exceed the growth velocities and circulation 
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patterns form in the melt. 

The approximate analytical solution can also provide considerable insight into 
the roles of the various processing parameters from which some general 
conclusions may be drawn: 

1. Axial segregation is diffusion controlled and radial segregation 
is directly proportional to the applied acceleration as long as the 
maximum buoyancy-driven flows are less than the growth velocity. 

2. The mixing between the radial flow and the diffusion boundary 
layer is the greatest when the diffusion length and the radial flow 
boundary layer thickness are comparable which results in the 
maximum radial segregation. The radial segregation becomes 
vanishingly small as either parameter becomes very small or very 
large. 

3. If the diffusion length is much larger than the ampoule radius, 
the radial segregation will be proportional to the growth velocity. If 
the ampoule is much larger than the diffusion length, radial 
segregation will be lessened by increasing growth velocity. 

4. If the diffusion length is much larger than the ampoule radius, 
the radial segregation will be a strong function of the ampoule 
radius, both because of the dependence of flows on the ampoule 
radius as well as the influence of sample size on mass transport. If 
the ampoule radius is much larger than the diffusion length, radial 
segregation will be influenced by sample radius only through the 
increased convective velocity. 

5. If the ampoule radius is much smaller than the diffusion length, 
the percent radial segregation is proportional to (1-k). If the 
ampoule is much larger than the diffusion length, the percent 
radial segregation is proportional to (1-k)/k. 

6. For periodic accelerations, the velocity amplitude is inversely 
proportional to the frequency for frequencies larger than 2.5 v/a 2 
Hz. The perturbation in composition at Hz the interface becomes 
proportional to the inverse square of the frequency for frequencies 
larger than 0.2 D/a 2 Hz, and proportional to the inverse cube of the 

frequency for frequencies larger than 2.5 v/a 2 Hz (assuming 
Schmidt number >1). 

7. Since the approximations for the governing equations used in 
this analysis were all all linear equations, any linear combination 
of solutions is also a solution. Therefore, the response to multiple 
frequencies, or to a combination of a steady acceleration and a 
number of applied frequencies can be found by simply adding the 
solutions for each individual case. 
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Although some of these conclusions may be self-evident or could have been 
drawn from the number of numerical computations that have been reported in 
the literature, it is useful and somewhat more assuring to be able to examine 
directly the interplay of the various parameters in the functional relationships 
that govern the process. 
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SECTION 3 


SOLUTE REDISTRIBUTION IN BRIDGMAN GROWTH 
WITH TRANSVERSE ACCELERATIONS 


3.1 INTRODUCTION 

The problem of obtaining compositional homogeneity in bulk crystals grown from 
the melt has received considerable attention in recent years. For example, 

Chang and Brown showed that even in vertical Bridgman growth with a stabilizing 
thermal gradient, buoyancy-driven convective flows resulting from radial 
temperature gradients will cause considerable radial segregation [1]. Application 
of strong magnetic fields are somewhat, but not completely, effective in 
controlling these flows in Earth's gravity [2], The problem is even more difficult 
when solutal gradients oppose the thermal gradient as was shown by Coriell, et 
al [3]. A number of attempts have been made to reduce these flows by taking 
advantage of various space flight opportunities to grow crystals in a microgravity 
environment. Witt and Gatos were able to eliminate dopant striations in their 
early experiments, but did not achieve uniform doping in a dilute system [4-5] . 
Fripp and Crouch observed almost complete mixing in their attempt to grow a 
non-dilute alloy-type system [6]. 

As a result of these early experiments, investigators began to realize that simply 
reducing gravity by 5 to 6 orders of magnitude may not be sufficient to eliminate 
solute redistribution due to residual convective flows. Chang and Brown's 
computations did indicate that 5 to 6 orders of magnitude reduction in gravity 
should be sufficient to virtually eliminate radial segregation in dilute systems up to 
approximately 1 cm in diameter, but their model only addressed steady 
longitudinal accelerations (g-vector along the furnace axis). More recently, 
Alexander, et al developed numerical models to treat transverse accelerations [7] 
and time dependent accelerations [8]. As might be expected, they showed that 
significant segregation could occur in dilute systems for transverse accelerations 
less than 1 micro-g. Because of this realization, an attempt was made to fly the 
first United States Microgravity Laboratory (USML-1) Spacelab mission so that 
the axis of the Crystal Growth Furnace (CGF) was aligned with the residual 
acceleration vector to minimize transverse accelerations. 

All of the models cited above utilize computational fluid dynamic (CFD) 
techniques to solve the coupled set of partial differential equations that describe 
the flow and compositional fields. While much can be learned from these 
computations, they provide data only at selected points in the multi-dimension 
parameter space involved in the process. Thus it is difficult to gain insight into 
how the systems would behave under different processing conditions or to 
extrapolate the results to other systems. Several previous attempts at 
developing scaling laws [9-10] using a single length scale to characterize the 
process were not particularly successful because of the complicated interplay of 
the various length scales that have to be considered. This was demonstrated by 
Naumann and Baugher who succeeded in developing an approximate analytical 
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model for dilute systems under steady and time-varying longitudinal accelerations 
which gave some insight into the complexities involved in scaling even a fairly 
simple system [11]. Garandet [12] used an impirical model similar to that used by 
Batchelor [13] to describe the flow in a 2-dimensional slot. He then used an 
order-of-magnitude scaling analysis to determine the lateral compositional 
variations as a function of the Peclet number, which he defines as the ratio of slot 
width to the diffusion length. Having obtained the functional dependence in the 
limits of high and low Peclet numbers, he used computational results to refine his 
order-of magnitude estimates and obtained an interpolation equation that gives 
quite good agreement with the computational results over the entire range of 
Peclet numbers. 

In this section an approximate analytical model is developed to describe the flow 
field and its effect on the concentration field in a Bridgman configuration 
subjected to very small steady transverse accelerations so that departures from 
diffusion limited behavior can be studied. Unlike the model developed by 
Garandet, which assumes a constant temperature gradient along the axis of the 
ampoule, a more realistic exponential density gradient will be assumed. This will 
allow non-dilute systems whose solutal fields decay exponentially along the axis 
to be considered also. The resulting equations also provide valuable insight into 
how the various parameters interact and should provide useful scaling laws that 
can be applied in general to different growth systems. Such scaling laws will 
serve as useful guides for developing future experiments and will allow more 
efficient use to be made of the numerical models. 


3.2. ASSUMPTIONS 

To keep the mathematics tractable, several simplifying assumptions must be 
made that limit the applicability of the model. First, the flow equations are solved 
independently of the mass and heat transport equations on the assumption that 
the flows are small enough that they do not significantly alter the density field. 
This is a good assumption for dilute systems where the density variations are 
primarily due to temperature differences since most materials of interest have low 
Prandtl numbers and heat transfer is predominantly by conduction for the flow 
regime of interest. For non-dilute systems, this assumption restricts the flow 
model to small perturbations in the concentration field. However, since the intent 
of the model is to examine the influence of convective flows near the diffusion- 
controlled regime, this is not a serious restriction. Having already limited the flow 
model to very small flow velocities, the non-linear inertial terms are second order 
and can be ignored which greatly simplifies the calculations. 

Second, no longitudinal gravity field is assumed, thus eliminating any effects of a 
stabilizing or destabilizing density gradient. This also eliminates any effect from 
radial temperature gradients which had been studied previously [1 1]. Again, for 
the dilute case, in which it has already been assumed that the temperature field 
is unperturbed by the small flows, this is a good assumption. In fact, Alexander, 
et al investigated several cases in which the acceleration vector was applied 
downward at 45° to the furnace axis [7]. For small accelerations, where the 
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segregation is directly proportional to acceleration, these cases produced the 
same result as that of a pure lateral acceleration with the same horizontal 
component. 

To avoid the complexities involved with non-axisym metric cylindrical coordinates, 
the flow and mass transport analyses are carried out in a two-dimensional 
rectangular slot. Alexander, et al used this simplifying geometry for the bulk of 
their studies but also carried out a full three-dimensional computation for a few 
cases for comparison. They found that the two models generally agreed to within 
20-25% of each other. Given the other simplifying assumptions that were 
necessary to make as well as the uncertainties in some of the thermophysical 
properties that go into such calculations, this accuracy appears adequate for the 
intended purposes of this model. 

Finally, it is assumed that the crystal is growing with a flat interface at constant 
velocity under steady-state growth conditions; i. e. the diffusion field is fully 
developed and does not change with time. The density field is assumed to be 
one-dimensional and is represented by an exponential function that decays with 
some characteristic length scale. The growth ampoule is assumed to be long 
enough so that neither the flow field nor the concentration field is affected by the 
end away from the growth interface. 

3.3 FLOW MODEL 


In order to estimate the effects of small transverse accelerations on the 
redistribution of solute in directionally solidified melts, a model for the flow field 
must first be established. There have been a number of attempts to obtain 
analytical solutions to describe the two-dimensional flow fields in horizontal 
chambers with differentially heated ends walls. [13-18]. For large length to width 
ratios, all of the models describe the horizontal flow in the core region as a one- 
dimensional cubic function symmetric about the centerline. This core flow is then 
matched to a more complex flow in the end regions. Unfortunately, the resulting 
equations for the flow in the end regions are generally too complicated to be 
useful for this type of analysis. The impirical polynomial representation used by 
Garandet assumes a linear density profile and therefore is not applicable to the 
problem at hand. 

The flow problem may be formulated by replacing u and v in the x- and y- 

momentum equations by u = - — and v = — to satisfy the continuity 

9y 9x 

equation, cross differentiating, and subtracting to eliminate the pressure terms. 
The result is a single 4 1h -order partial differential equation for the stream function 

vjr. Neglecting the inertial terms, this equation becomes 


. o -jV | a4 V | = q 

dx 4 3x 2 9y 2 9y 4 ji 9x p 9y 


(3.1) 


104 



where g* and cfy are the components of the gravity vector and p(x,y) is the density 
of the fluid. It is assumed that the system is aligned such that the acceleration is 
predominately in the -y direction (horizontal Bridgman configuration) or that the 
density variations in the y direction are small enough that 


9y 0p 
p 3x 


» 


5x_ ^P 

p 3y 


(3.2) 


which could apply to a dilute system in a slightly misaligned vertical Bridgman 
configuration provided the thermal conductivity of the melt is high enough so that 
the isotherms are not affected by the resulting convective flows. 


The boundary conditions require that \p and its normal derivative vanish at each 
wall to enforce the no-slip condition. 

An exponential function of the following form is chosen to represent the density, 


p(x) = 


(P(0) ~ P(L))e~ xa - p(0)e~ L/x + p(L) 
1 - e' LA 


(3.3) 


where x = 0 is taken at the solidification interface and X is the e-folding length. If 
L is finite, this function describes the density profile as a decaying exponential 
between p(0) and p(L) with a gradient given by 

I--*?- «» 

Introducing dimensionless coordinates £ = x/a, r| = y/a, equation (3.1) can be 
written. 


a 4 y o a> a> _ g[p(Q)-p( oo )]a 4 e~ a ^ _ 

a^ 4 a^an 2 an 4 p 0 vx 


(3.5) 


where Gr = 2 ^ anc * A P is 9i v © n by eq. (3.4). 


(Note that if L » X, Gr ^] a . whereas if L « X, Gr 

P(L)v X. 


g[p(0)-(L)1a 4 


P(*-)v 2 L 
Boundary conditions are: 


and p(x) becomes linear from 0 to L.) 
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y=0 at £ = 0,°° and rj = ± 1 ; 


v = — ^ = 0 at q = 0,oo; and u= — --^ = 0 at n = ±1. 

a 3ri 

These boundary conditions assume a semi-infinite ampoule. They could be 
enforced at some finite length, but at considerable complexity. However, since 

the region of interest is in the vicinity of £ = 0, and since the comparison with 
computer modeling is good provided L » a, the above formulation appears 
adequate for the intended purpose. 

If the gradient were uniform over a distance L»h, the equation would be 
dominated by the r| -derivatives in the central regions away from the ends and the 
stream function behaves as a simple 4th deg. polynomial in 1 ]. We might 

therefore assume that the stream function preserves this form for all \ and write 
an approximate solution in the form 

V = <P(0(t| 4 -2t| 2 + i), 0 < ^ < L/a (3.6) 


where tp(£) is a function of £, only and represents the stream function along tj = 

0. Putting this back into eq. (3.5) and evaluating the r\ derivatives at t] = 0, we 
obtain the 4th order ordinary differential equation 


cT 



+ 24 cp = Grve -a5 '\ 


(3.7) 


This may be solved by Laplace transforms using the boundary conditions <p(0) = 
<p'(0) = 0 and rp(^) = <p'(£,) -> 0 as £, -» «> to obtain 


<p(S) = Gr*v je a ^ /A ’-[cos(p^) + ysin(p^)]e 

where 


Gr’ = 


Gr 

(24-8a 2 /^ 2 +a 4 / A. 4 ) ’ 


a 2 = 2+V6 , p 2 =-2+V6, 


(3.8) 


and y = 


(a-a/X ) 

P 


The complete stream function is therefore, 


\| f = Gr‘v je a ^ a -[cos(p£)+ysin(|3£)]e }(ti 4 — 2 tj 2 + 1) . (3.9) 

The horizontal velocity may be obtained by differentiating 
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U = “af^ =_4Gr ta){ e a ^M cos (f^) + T sin (l^)]e °^}(ti 3 - r|) . (3.10) 

Likewise the vertical velocity can be found from v = — 

a 


v=Gr*| — )<!-— e -a ^ /A - 
a 


a 0 a ^ a + [(a-py)cos(p^) + (ay + p)sin(p^)]e < ^|(r| 4 -2 ti + l) .(3.1 1) 


3.4 COUPLING FLOW MODEL TO MASS TRANSPORT 

In a reference frame that moves along the x-axis at the growth velocity v g so that 
the growth interface remains stationary in the lab system, the concentration field 
C is given by the two-dimensional, steady state mass transport equation 


( 



ac dc 

dx dy 


= D 


a 2 c a 2 c 

dx 2 dy 2 


(3.12) 


where D is the diffusion coefficient. 


Assuming a plane front growth interface, the boundary conditions are 
BC^\ 

D ~faj + v g (1-k)C(0,y) = 0, C(°° y ) = C„ for ally 

and 

D^~) =0 for all x 

Jy=±a 

where k is the segregation coefficient. 

The concentration field C can be written as 


C(x,y) =C 00 (C x (x) + C(x,y)) (3.13) 

where C„ is the bulk concentration of the melt, is the solution to the steady 

state, purely diffusive mass transport equation, and C(x,y) is the correction from 
convective flows in the melt. 

An expression for the first order correction from the convective flows can be 

obtained by assuming that u, v, and 6 are small so that the second order terms 

(products of u or v and derivatives of C) can be ignored. Dividing the mass 
transport equation by v g , inserting eq. (3.13), and retaining only first order terms, 
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(3.14) 


sT a 2 c a 2 c ~ 

^ a a^ 2 an 2 

where 8= D/v g and u = u(x,y)/v g . 

The Cx(x) is given by 

C x (x) = (1/k-1)exp(-x/5)+1. 

Let 

U = U 5 (« ^(Tl) 


(3-15) 

(3.16) 


and assume that C(x,ri) can be written as 


C(5,t l) = f(5)exp(-a£/S )g(ri) 


(3.17) 


where f(^) and g(r]) are respectively functions of ^ and q only. 

Taking derivatives of eq.(3.15) and eq. (3.17), putting them back into eq.(3.14) 
and rearranging terms, the equation can be written 



a df(5) 
S 3E, 


(1/K-DSlMM 

v ' 6 2 g(Ti) 


f (5) a 2 gh) 

g(Ti) an 2 


= o. 


(3.18) 


If the (rj) had the same functional form as g(r|), the variables would separate 
and g(q) could be expressed as a Fourier series. However, the boundary 

conditions at q = +1 require 3g/dr| = 0 whereas 9u / 9ri * 0. Thus it is not 
possible to obtain a mathematically consistent solution by separation of variables. 
Guided by the numerical results, it would seem that the ^-dependence of C could 
be reasonably represented by g(r|) = -sin(7cn/2). This satisfies the boundary 
conditions at t| = ±1 and, like the (q 3 - ti) term in the expression for u, is 

antisymmetric with respect to r| = 0. To obtain an approximate solution, the u 
term in the mass transport equation will be represented by 


u = -|u 5 (4)si n (T|n / 2) 


(3.19) 


in which the (r| 3 - r\) term has been replaced by -jc/8 sin(nji/2). The coefficient jt/8 
was chosen to give the same integrated flow as before, i.e. 
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1 1 

J ("n 3 — Ti)cfri = - J|sin(T|7i/2)dr|. 

0 

Since this flow acts only on the gradient of Cx, which is not a function of r), the 
artificial distortion introduced into the flow field for mathematical expediency 
should not produce a serious error. 

The 2nd order partial differential equation can now be reduced to an ordinary 
second order differential equation which using eq. (3.10) can be written 

f ( 5 ) “ ( 5 ) f(5) = G je -a ^ a - [cos (p ^) + y sin((3 ^)]e _a ^ j (3.20) 

where G = (l/k-l)^ J Gr v/a = (1/k -1)— Gr*Sc— 

5 2 2 v g v '2 5 

and Sc is the Schmidt number (v/D) 

The solution has the form 


f(0 = G{A 1 e" a ^ + [A 2 sin(K)+A 3 cos(K)]e-* +A 4 e^} 


where 


A 1 = 


- 1 


a 2 a 2 jt 2 

xF + wPT 


a 2 


yE-f a 

E 2 + F 2 ’ 3 


yF+E 
E 2 + F 2 ’ 


(3.21) 


E = 4+ aa/S- 7t 2 /4, F= (3(2a +a/5) . 


The quantity r is the negative root of the indicial equation, r 2 -ra/8-rc 2 /4 = 0 
which is given by 


r =^( 1 -VT+?5 2 7a 2 J. 


(3.22) 


(Since the perturbation must die out at large x, the second unspecified coefficient 
corresponding to the positive root of the indicial equation is set to zero.) 

The remaining free coefficient, A 4 , must be determined from the conservation of 
solute boundary condition at the interface, i.e., 


D 


dC^ 

dx Jx = o 


+ v g (1-k)C(0,y) = 0. 


Inserting eq.(13) and eq.(17), this condition requires 


(3.23) 


109 



3.5 COMPARISON WITH NUMERICAL COMPUTATIONS 

The approximate analytical model was extensively tested against numerical 
computations. The results shown are from the following test case. 
g = 980x1 O' 6 cm/sec 2 

pAT= 0.025 
D = 0.00013 cm 2 /sec 
v = 0.0013 cm 2 /sec 
L = 2.5 cm 
a = 0.5 cm 
X = 1 .0 cm 
k = 0.1 

vg = 6.5 (im/sec 

The stream function was computed numerically from eq. (3.5) using an equally 
spaced 21x51 grid and a simple relaxation method. For this computation, the 

boundary conditions \|/(r|,L/a) = \y’(T|,Lya) = 0 were imposed instead of allowing 
these functions to approach 0 as x goes to °° as was done for convenience in the 
analytical model. Comparison of the numerical results against the stream function 
calculated from eq. (3.9) is quite good as may be seen from Figure 3.1 which 
compares the stream function along y = 0 and from Fig. 3.2 which compares the 
stream lines. The amplitude and location of the peak agree to within a few percent, 
the slight difference apparently results from the different boundary condition 
imposed at x = L. The most critical aspect of the model is its ability to predict the 
horizontal velocity profile in the region of the interface. This it does extremely well 
as may be seen in Fig. 3.3 which compare the horizontal velocity profiles 
calculated from eq.(3.10) with the numerical result. Fig. 3.4 compares the 
maximum horizontal and vertical velicity profiles. Similar comparisons were made 

for different values of X as shown in Fig. 3.5 and 3.6. As can be seen, the model 

improves with increasing X. The agreement remains quite good for values of sJX < 

1 where the viscous terms in the core region are dominant, but the model starts to 
break down for a IX > 1 because the assumption that the y-dependence of the 
stream function in the core region holds everywhere is apparently not valid near 
the corners. However, the model is still within a factor of 2 of the numerical result 
as a/X get large. 

The composition profile in the melt was also computed numerically using the same 
relaxation method. The isoconcentrates calculated by combining eq. (3.15) and 
eq. (3.26) using the coefficients defined in eq.(3.21) are compared with those 
computed numerically in Fig. 3.7. The y-dependence of the concentration at the 
interface is shown in Fig. 3.8. Again the results are virtually identical. 

Having demonstrated that the analytical model can faithfully reproduce the 
essential effects of a small lateral acceleration, comparisons with the results of 
Alexander, Ouazzani, and Rosenberger (AOR) [7] were then attempted. A choice 
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dm 

% 


(3.24) 


■1 -^f(o) = o. 

k=o 6 

A 4 may now be found using eq. ( 21 ) 

A = (ka/S+a/X,)A-i 1 - pA 2 +(ka/5+oc)A 3 
4 r-ka/5 


(3.25) 


Finally, the perturbation to the concentration field is obtained by combining 
eq.(3.17), (3.19) and (3.21) to obtain 

= G e_a? /s { + [ A 2 sin(p^) + A 3 cos( p£)]e““^ + A 4 e r ^ } sin(r| 71 / 2 ). (3.26) 

The concentration at the interface is then 

~ = 1/k + G (A 1 + A 3 + A 4 )sin(ri7i/2). (3.27) 

The segregation along the growth interface is characterized by a segregation 
parameter £ defined as the difference between the maximum and the minimum 
concentration divided by the average composition along £= 0. Since the maximum 
and minimum concentration will occur at r| =+1 and -1 respectively, and the 
average concentration will be CJW, we can write 




C(o, + i)-C(o.-i) 


C / k 


=|2kG(A 1 + A 3 + A 4 )| 


Using the definition for G in eq. (20), this may be written as 




, GrSc(a/8)(A 1 + A 3 + A 4 ) 

' ' r\ a o 1 \2 . / _ \4 


24-8(a/X) +(a/X.) 

The quantity At + A 3 + A 4 can also be written as 


A + A + A = + PA 2 +(r+a)A 3 

1 3 4 r-ka/5 


(3.28) 


(3.29) 


(3.30) 
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Stream function along y=0 



Dimensionless length (x/a) 

Figure 3.1 . Comparison of the approximate solution for stream function along the 
x-axis with a numerical solution to the flow equations for the test case. The 
numerical computation was carried out for a finite (2.5 cm long) ampoule, whereas 
the approximate solution assumed an infinite ampoule. Note the magnitude and 
shape of the stream function are accurately captured by the analytical model; the 
major differences being in the boundary conditions at the end. 
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model. Contours are 2 x 1 0-® cm 2 /sec. 




Scaled width (y/a) 


Scaled width (y/a) 


Stream Function (Numerical) 




Figure 3.3. Comparison of horizontal velocity field for the test case computed numerically (top) with the approximate analytical model. 
Contours are 1 x 10 -5 crr^/sec. 



Horizontal velocity (numerical) 




Umax and Vmax 



Dimensionless length (x/a) 


Figure 3.4. Comparison of the maximum horizontal and vertical velocities 
(dashed lines) calculated from the approximate analytical model with those 
computed numerically (solid lines). 
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Figure 3.5. Comparison of stream function along t\ = 0 computed numerically 

(solid curve) with the approximate model (dashed curve) for X = 0.25 cm, L= 2.5 
cm, a = 0.5 cm. 
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0.006 



Length ( cm) 

Figure 3.6. Comparison of stream function along r\ = 0 computed numerically 
(solid curve) with the approximate model (dashed curve) for X = 10 cm, L= 10 
cm, a = 1 cm , Gr v = 1 .582 cm 2 /sec. 
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Figure 3.7. Comparison of concentration profiles in the melt for the test case computed numerically (top) with the approximate analytical model. 
Contours are 2, 3 10 C„ 



Scaled Width (y/a) Scale width <y /a > 


Concentration Profiles (numerical) 




Melt Composition at Interface 



Figure 3.8. Comparison of the compositional profiles in melt at the growth 
interface for the test case. The solid curve was computed numerically and the 
dashed curve was calculated from the approximate analytical model. 
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of a JX = 2 produced a density profile that was reasonably similar to that produced 
by the average axial temperature profile used in the AOR computations. As it turns 
out, the segregation is not particularly sensitive to the choice of aA. in this regime 
as will be shown later. 

As may be seen in Table 3.1 , the agreement is quite good as long as the 
perturbation is small (i.e. £ < 1) which corresponds to Gr Sc less than -100. Within 
this range the conditions required for the analytical model to be valid are generally 
met; i.e., the flows do not perturb the density field significantly and the second 
order terms in the mass transport equation are small enough to be safely ignored. 
(In fact, the agreement is quite good beyond this range, but since the perturbation 
is no longer small, this agreement must be considered fortuitous.) Hence, for Gr Sc 
< 100, the degree of segregation is directly proportional to the product of Gr Sc (for 

fixed a/5). At first glance, this appears to contradict the findings of AOR who argue 
that the segregation does not depend linearly on Sc. It should be remembered that 
this model is restricted to small enough flows so that nonlinear effects have not yet 
manifested themselves. Also, AOR altered the Sc by varying the diffusion 

coefficient, keeping the growth rate constant. This of course would cause 5 to 
change, thus altering the coefficients defined in eq. (3.21). Actually, for constant 
values of a/8, the AOR results appear to scale linearly with Gr Sc up to -360. 
Beyond this point the convective mixing becomes such that the segregation is 
reduced. 


Table 3.1 


CO 

O 

a 

(cm) 

m 

AT 

Gr Sc 

a/5 



10‘ 4 

0.50 

6.50 

100 

3600 

2.50 

0.80 

9.54 

IO - 5 

0.50 

6.50 

100 

360 

2.50 

0.927 

0.954 

10' 5 

0.50 

0.65 

100 

360 

0.25 

0.119 

0.181 

10-5 

0.25 

6.50 

100 

45 

1.25 

0.120 

0.089 

5x10-5 

0.50 

6.50 

100 

1800 

2.50 

0.542 

4.77 

10' 6 

0.50 

6.50 

100 

36 

2.50 

0.113 

0.095 

10‘ 6 

0.50 

0.65 

100 

36 

0.25 

0.02 

0.018 

1.41x10-5 

0.50 

6.50 

115 

590 

2.50 

1.52 

1.56 

1.41x10-6 

0.50 

6.50 

115 

59 

2.50 

0.215 

0.156 

1.41x10-6 

0.50 

0.65 

115 

59 

0.25 

0.015 

0.024 

10-5 

0.50 

6.50 

20 

72 

2.50 

0.226 

0.191 

10-5 

1.00 

6.50 

20 

580 

5.00 

0.645 

1.156 

IQ' 6 

0.50 

6.50 

20 

7.2 

2.50 

0.023 

0.019 
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3.6 COMPARISON WITH GARANDET’S RESULTS 

Garandet considered the case of a linear thermal gradient which can be duplicated 
in this model by letting X For small Peclet numbers, his order of magnitude 
analysis obtained a scaling law that estimates the variation in concentration across 
the interface AC to be given by 

AC=Pe | G ^ c(1 ~ k) C„, Pe'«1 (3.31) 

768 k 


where his Peclet number Pe’ = v g H ID and Grashof number Gr’ = (3j g (AT / L) H 4 / 
v 2 . Comparing this with computer analysis, he adjusted this scaling law by 1/2.3 to 

obtain agreement with the numerical results. Since £ is defined as k AC /C „ , 
Garandet’s adjusted result can be written in terms of the Gr and Pe defined with 
length scale a instead of H as 


C = 


kAC _ GrSc(1-k)(H/a) 5 
C„ “ (2.3)768 


Pe = 0.0181(1-k)Gr Sc Pe 


(3.32) 


For X — » oo, the segregation parameter given by eq.(29) and eq.(30) reduces to 

(3.33) 


r . .GrSc_ 

£ = jc (1 — k) Pe 


24 


rA 1 -(3A 2 + (r+a)A 3 
r-kPe 


where Pe = a/8. 


The coefficients Ai, A2, and A3 also simplify and can be written as polynomials of 
Pe. For small Pe, only the constant terms need be considered. Inserting the 

numerical values for a = ^2 + V6 and (3 = ^/-2+V6 , these coefficients are given 
by 


A, 

A 2 

a 3 


4 / 7t 2 =0. 405 

(g/p)(4-»*/4)-2gp 018 

(4-ti 2 /4) 2 + 4a 2 (3 2 

2 a 2 + 4 r ,t2/4 = 1.008 . 

(4-7i 2 /4) +4a 2 p 2 


Also for small values of Pe, the r in eq. (3.22) = - tc/2 which is » k Pe. Inserting 
these values into eq. (3.29) and (3.30), we obtain 

£ = 0.019 (1 - k) Gr Sc Pe , Pe «1 (3.34) 
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which is virtually identical with Garandet’s result. 

For large Pe, Garandet obtains a scaling law from his order of magnitude analysis 
that predicts 


AC = 


1 Gr'Sc (1-k) 
Pe' 3 48 k 2 


Pe'»1. . 


(3.35) 


He finds that he must adjust this result by a factor of 12 to bring it into agreement 
with numerical results. With this adjustment, the segregation parameter predicted 
by his model can be written as 


kAC _ (1-k)GrSc(H/a) 1 
C m ~ k 48 Pe 3 


0.500^— ^GrSc-^ 
k Pe 3 


(3.36) 


Expanding eq. (3.30) in terms of 1/Pe and inserting the numerical values for a and 
P, the 1/Pe and 1/Pe 2 terms vanish leaving 


A^ ^3 ^4 


r A, ~PA 2 +(r+g) A 3 ^.BBSQ/Pe^ + Otl/Pe 4 ) 

r-kPe -kPe 4 ' ' 


Putting this into eq. (3.29), the segregation is given by 


y / h , \ GrSc _ 

C = 7i(1-k) Pe 

s v ' 24 


-4.337(1/Pe 3 ) + 0(1/Pe 4 ) 
-kPe 


= 0.639^— ^Gr Sc -L .(3.38) 
k Pe 3 


It may be seen that the present model not only gives the same scaling as obtained 
by Garandet for the case of X -» °° , but also gives very close agreement with his 
numerical results without any adjustments. This attests to the accuracy of the 
model. 

3.7 EFFECT OF THE EXPONENTIAL GRADIENT 

Since neither this model nor Garandet’s model imposes boundary conditions on 
the fluid at the end of the ampoule away from the growth interface, the length L can 
represent a linear density profile in a semi-infinite ampoule. Figure 3.9 plots the 
ratio of the segregation resulting from an exponential density profile to that of a 
linear profile where the densities at 0 and L are the same. (L in this case was 

taken to be 5 a.) Increasing values of a \TK represent the steepening of the initial 
density gradient. The effect of the exponential density profile is to enhance the 
degree of segregation as the initial gradient steepens up to A. « a. The effect is 

more pronounced for larger values of Pe (smaller 6) because more of the diffusion 
field is in the steeper part of the gradient where the convective flow is stronger. As 
the initial gradient steepens beyond a certain point, however, the flow becomes 
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2.5 



a/lambda 

Fig. 3.9 Ratio of the segregation resulting from an exponential density profile to 
that resulting from a linear density over a distance L = 5 a for k = 0.1 . The 

densities at 0 and L were taken to be the same. alX =1 represents a linear 

gradient which steepens as a/X increases. The Pe = v g a /D or the ratio of 
ampoule width to diffusion length. 
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more confined to the region of the solidification interface which causes the velocity 
to diminish because of the higher shear forces as may be seen in Fig. 3.10. 


The interplay of the steepness of the gradient and the Peclet number can be seen 
in Figure 3.1 1 . The contours represent the segregation £ scaled by the product Gr 
Sc (A/a). This effectively scales the segregation by a form of the Gr that does not 
involve A (i.e. g (Ap/p) a 3 / v 2 ) so that now as a/A 0, the gradient vanishes. It 
can be seen that the segregation is maximized for Pe ~ 2 for all values of a/A. The 
segregation increases with increasing a/A up to a/A = 1 .4. 

For non-dilute systems where the convective flows are driven primarily by the 
solutal gradient, A = 8 and a/A = Pe. The interplay of Pe and k for this situation is 
shown in Fig. 3.12. Note that the Pe that produces the maximum segregation 
decreases from 2 for k = 0.1 to 0.8 for k = 10. The interplay of the streamlines and 
concentration gradient for these two cases under conditions of maximum coupling 
may be seen by comparing Fig. 3.13 and 3.14. A similar result is obtained for the 

non-dilute case in which a/A is fixed. 


3.8 SUMMARY AND CONCLUSIONS 

In order to study effects of small sustained transverse accelerations on the flow 
and concentration fields in a Bridgman-Stockbarger crystal growth apparatus, the 
molten part of the system has been modeled as a fluid in a semi-infinite two- 
dimensional horizontal slot with a density field that falls off exponentially from one 
end. It is assumed that the flows are sufficiently weak that the density field is not 
affected by the flow, thus allowing the flow equation to be solved independently of 
the mass transport equation. The maximum flow velocity was found to be ~ 

Gr (v/a)/9V3 . Therefore, for dilute systems, in order for the thermal Peclet 
number (v A / k) < 1 , the thermal Rayleigh number g (Ap/p) a 3 / v k < 1 5. Similarly, 

for non-dilute systems, the solutal Rayleigh number g (Ap/p) a 3 / v D < 15. 
Stabilizing or destabilizing effects of a longitudinal acceleration component have 
not been considered in this work, although such effects could be extremely 
important, especially for non-dilute systems. 

A closed form approximate solution was obtained for the two-dimensional flow 
equation by assuming that the cubic y-dependence of the velocity field in the core 
region holds for all x. This allows the biharmonic partial differential equation for the 
stream function to be reduced to a fourth-order ordinary differential equation which 
was solved by the method of Laplace transforms. The resulting solution is a 
relatively simple expression which can be differentiated to obtain the horizontal and 
vertical velocity fields. Although this is only an approximate solution, excellent 
agreement with numerical solutions was obtained provided the gradient region 
was longer than the half-width of the chamber. This covers most of the region of 
interest for this application. Even for gradient lengths less than the chamber half- 
width, the model was within a factor of 2 of the numerical result. 
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Figure 3.10. Maximum horizontal velocity profile for X = 0.1 a (dotted line), X = 
0.5 a (dashed line), X = 0.8 a (solid line), X = 2 a (dash-dash-dot line), and X = 1 0 
(dot-dash line). The velocity is scaled by Gr (v/a). 
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Fig. 3.1 1 . Contour plot of segregation scaled by GrSc(A/a) for k = 0.1 . Peak 
value is 0.00566 at afk = 1 .25 and log Pe = 0.3. Contours are spaced at intervals 
of 0.0005. 
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- 2-10 1 2 

Log Pe 

Fig. 3.12. Segregation scaled Gr Sc (A/a) for the case of solutal-driven convection 
in which a JX = Pe. Peak value is 0.096 and contour spacing is 0.001. Note the 
shift in the Pe that causes the maximum segregation as k goes from 0.1 to 10. 
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Figure 3.13. Concentration field overlaid on stream function for maximum 
coupling (Pe = 2.0) between the flow and the concentration gradient for k = 0.1 . 
Isocontours of the concentration field represent 2,3.... 13 C^. 
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Figure 3.14. Streamlines and Isoconcentrates for maximum coupling (Pe = 0.8) 
between flow and solute field for non-dilute case with k = 10. Contours are (from 
the right hand side) 0.9,0.8,....C eo . 


The resulting expression for the density-driven convective velocity field was then 
substituted into the mass transport equation and a first-order perturbation solution 
was obtained for the solute redistribution. This solution is valid so long as the 
second order terms (i.e. products of the convective velocity and gradients of the 
perturbed solute concentration field) can be ignored. Since first order terms 
involve products of the convective velocity and gradients of the unperturbed 
density field, this condition is met so long as the perturbation is small. This 
restriction is not serious since the primary purpose of the study was to examine the 
departures from diffusion-limited concentration fields caused by small transverse 
accelerations. Again the results from the perturbation solution were found to 
compare favorably with numerical computations. 

With this analytical model, the radial segregation produced by small transverse 
accelerations could be analyzed as a function of chamber width a, the diffusion 

length 8, and the e-folding length of the exponential gradient X. The model was 
compared with an earlier model developed by Garandet who assumed a linear 
density profile. In both models the degree of segregation scales as the product of 

Gr Sc Pe (1 - k) for Pe =a / 5 «1 and as Gr Sc (1 - k)/ k Pe 3 for Pe »1 . By setting 
X to oo to obtain a linear density profile, the present model produced virtually the 
same results obtained by Garandet after he had adjusted his scaling analysis to 
agree with numerical computations. The effect of steepening the gradient 

(shortening X) while keeping the density the same at x = 0 and at x = L is to 
increase the degree of segregation by up to several fold, depending on the Pe. 

For k = 0.1, the segregation parameter (maximum composition variation along the 
width of the sample divided by the average composition) can be as high as 
0.00567 Gr Sc (X/a) under conditions of maximum coupling between the flow field 
and the density field. If compositional variations are to be held to less than 1 %, the 
product Gr Sc (X/a) < 1 .76. Since Sc numbers are generally 0(10) for 
semiconductors, Gr (X/a) = g (Ap/p) a 3 / v 2 < 0.176. For a typical system such as 

the Ga-doped Ge modeled by Alexander et al, (Ap/p = 0.025, v = 0.0013 cm 2 /sec, 
a = 0.5 cm), the transverse acceleration must be less than 0.1 micro-g for 
conditions of maximum coupling. This illustrates the extreme sensitivity of such 
systems to sustained lateral accelerations. 

The situation can be even more critical for non-dilute systems because of the 
greater density differences due to solutal gradients. Furthermore, many of the ll-VI 
systems of interest (e.g. Hg x -iCd x Te, Hg x .iZn x Te, etc. have Sc numbers O(10 2 - 
10 3 ). To achieve compositional uniformity in such systems by simply reducing the 
effective gravitational acceleration without regard to direction of the steady or 
quasi-steady component of the residual acceleration, one would have to have 
residual accelerations considerably less than 0.1 micro-g to assure diffusion limited 
solute redistribution in any growth system of reasonable scale ( = 1 cm in diameter). 
By understanding the interplay of the various length scales, some relief can be 
obtained by de-tuning the experiment away from the most sensitive set of 
conditions. However, in a non-dilute system of given diameter, the only adjustable 
parameter is the growth velocity. Here, the experimenter has a limited range in 
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which to work since he must consider the time available to grow a usable amount 
ot material if he wishes to grow at a slower rate, while faster growth would require 
a higher thermal gradient to stabilize the solidification interface which may not be 
desirable or achievable. 

Therefore, it is clear that simply reducing the effective gravitation acceleration by 5- 
6 orders of magnitude cannot be expected to achieve diffusion limited growth 
conditions for the systems of interest. Extreme care must be exercised in 
controlling transverse accelerations by orienting the furnace axis along the residual 
acceleration vector. The possible benefits of maintaining a small continuous 
stabilizing acceleration should be explored as well as the application of a modest 
axial magnetic field to further damp the residual flows. 
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SECTION 4 


PREDICTIONS OF EFFECTS OF RESIDUAL ACCELERATIONS ON USML-1 

The effects of residual accelerations on Bridgman growth of dilute systems can 
be estimated from the following: 

For accelerations along the axis of the furnace, the maximum convective velocity 
is given by [1], 


w = - 


g y(at) i 

48 v 


(4.1) 


The perturbation in the concentration field at the growth interface is 


C(r,0) = PeX E„J 0 (rs„)f„(0) 

n=1 


(4.2) 


where the Peclet number Pe=w/ v_, 


En— " 


64 


Ji(as n ) L(as n ) 3 (as n f 


and 


2 G 

p 3 + r n 

3 p 2 + 3 a p + a 2 + s n 2 , 

+ Sn 2 (3 p + a ) 

r n - ka 

. 

,p 2 +a p - s n 2 , 



The parameters a, (3, and s n are reciprocal lengths given by 


oc = D/Vg , p — VX, and Sn — (4.3) 

where £, n are the zeros of Jo and X is the thickness of the momentum boundary 
layer of the flow across the face of the growth interface. This is estimated to be 
half the distance to the point where the density difference is balanced. For 
thermally driven convection, it is assumed that the thermal profile is given by 

T(z) = (T.-T m .,)(1-e- z/5 )-T„, (4.4) 

where the thermal length 5 = (T„ -T me „)/the imposed thermal gradient. The X 
is then given by X = 8 ln(2)/2. 

The quantity r n is given by 
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V 1 +4 s n 2 /a 2 ! 


(4.5) 


r„=Ml- 

2 


and G is 


G--(l-k)a 2 p 2 ei— . (4.6) 

4 k 

The effects of lateral acceleration may be found using a modification of the 
procedure outlined in [1]. In this case the maximum velocity is given by 


.. = y(T ho , T^Ja g ( 4 - 

max 9V3(a + 2^)v ' 

The maximum perturbation in the concentration field at the growth interface is 
given by 


C(r,0 ) = Pe £ F n sin(2 n n r/a) f n (0) 

n=1 


(4.8) 


where 


F n = 


-18V3(-1) n 


7t 3 a 3 


4.1 Se-doped GaAs 

For the growth of Se-doped GaAs (Dave Matthiesen’s experiment) the following 
are assumed: 


a = 0.75 cm 
y= 5x10- 4 /K 
T hot = 1255°C 
Tmelt = 1238°C 
grad T) 0 = 10 K/cm 

<AT> = 1.7 K 
k = 0.1 

v = 2.98x1 0'3 cm 2 /sec 

D = 1x1 O' 4 cm 2 /sec 

v g = 2.5x1 0‘ 4 , 5x1 0‘ 4 cm/sec. 
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From the thermal profile, the momentum boundary layer of the radial flow was 
estimated to be X = 0.589 cm. The computed maximum convective velocity is 3.3 
(g/go) cm/sec for acceleration along the furnace axis. For v g = 2.5x1 O' 4 cm/sec, 

the maximum perturbation is C = 497(g/g 0 )C„ /k which implies the radial 

segregation will be C = 497 (g/go). Therefore, if £ < 1%, the maximum steady 
acceleration along the furnace axis must be < 20 micro-g’s. 

For the second growth rate, v g = 5x1 O' 4 cm/sec, £ = 577 (g/go)- Therefore, the 
maximum steady acceleration along the furnace axis must be < 17 micro- 

Q’s. 

For transverse accelerations, the u ma x = 39 (g/go) For v g = 2.5x1 O’ 4 cm/sec, the 

maximum perturbation amplitude is C=7000(g/g 0 )C„, /k. For this case the 
difference between the maximum and minimum composition is twice the 
perturbation amplitude. Thus the radial segregation will be £ = 1.4x10 4 (g/go)- 

Therefore, for £ < 1%, the transverse acceleration must be <0.7 micro-g. 

For the second growth rate, v g = 5x1 0"4 cm/sec, £ = 1.7 xIO 4 g/go. Therefore, 
the maximum steady transverse acceleration must be < 0.6 micro-g’s. 

Several of the physical properties that were needed in the calculation are not well 
known, e.g.. y and D, and had to be estimated. It is believed, based on values for 
similar systems, that these estimates are probably within a factor of 2 of the 
correct value.. 

This experiment has fairly modest requirements for acceleration along the axis 
because of the low thermal gradients in the melt. Even so, the radial distribution 
is extremely sensitive to transverse accelerations owing to the fairly large 
diameter and small distribution coefficient. Since the melt is confined by a spring 
loaded plunger, there are no free surfaces and, as was shown in [1], the effect of 
periodic accelerations falls as the inverse cube of the frequency for frequencies 
greater than 10‘ 4 Hz. Therefore, this experiment will be virtually immune to 
crew-induced vibrations on the Shuttle. 


4.2 CdZnTe 

For the Cd.gsZn 0 sTe experiment (Dave Larson’s experiment), the following 
parameters were assumed: 
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a = 0.75 cm 
7=5x1 0'4 /K 
T hot = 1 1 15°C 
T me tt = 1095°C 
grad T) 0 = 10 K/cm 

<AT> = 2.0 K 
k= 1.2 

v = 4.35x1 O' 3 cm 2 /sec 
D = 1 xl O' 4 cm 2 /sec 
v g = 4.4x1 O' 5 cm/sec 

The small mole fraction of the ZnTe component and the near unity distribution 
coefficient for this system allow it to be approximated as a dilute alloy in which 
buoyancy forces are dominated by thermal rather than solutal effects. From the 
thermal profile, the momentum boundary layer of the radial flow was estimated to 

be X = 0.695 cm. For steady accelerations along the axis, the computed 
maximum convective velocity is w=2.64(g/ g 0 )cm/ sec. The maximum 

perturbation is C = 17(g/g 0 )C„ /k which implies the radial segregation will be £ = 
17 (g/g 0 ). Therefore, if £ < 1%, the maximum steady acceleration along the 
furnace axis must be < 588 micro-g’s. 

For transverse steady accelerations, the maximum convective velocity is 
estimated to be 28 (g/go) cm/sec which gives a maximum perturbation in the 

concentration field at the growth interface C=208(g/ go)CL /k. The radial 

segregation is C, = 416 (g/go)- Therefore, for £ < 1%, the transverse 
acceleration must be <22 micro-g. 

The modest temperature gradients, low solute concentration, and near unity 
distribution coefficient (k value) for this experiment combine to require very 
modest acceleration requirements, especially along the furnace axis. There will 
be a small void in the melt to allow for thermal expansion will may make the 
system more susceptible to transient and periodic accelerations (g-jitter), but the 
relative insensitivity of the system to steady accelerations will serve to reduce 
such effects. Therefore, it is not expected that any effects from non-steady 
accelerations will be observed. 


4.3 HgZnTe 

This is a non-dilute system in which both thermal and solutal gradients drive 
significant convective flows. For the case of axial accelerations, the solutal 
gradient and the radial thermal gradient responsible for radial segregation are 
perpendicular to one another. This produces a coupling between the perturbed 
concentration field and thermally driven convection. Since the model [1 ] for 
radial segregation caused by accelerations along the furnace axis does not 
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consider this coupling, it is no longer valid. However, for transverse 
accelerations, the solutal and thermal gradient are in the same direction and for 
small flows the model can still be used with minor modifications. 

For the Hg 84 Zn 16 Te experiment (Alex Lehoczky’s experiment), the following 
parameters were assumed: 
a = 0.4 cm 

y=5x10- 4 /K 
T hot = 800°C 
T melt = 695°C 
grad T) 0 = 60 K/cm 

<AT> = 15 K 
k = 4.35 

v = 8.0 xl 0'3 cm 2 /sec 
D = 6x1 0" 6 cm 2 /sec 
v g = 4.1x1 O' 6 cm/sec 
PHgTe = 8.0 
PZnTe = 5.35 


The formula for u max must be modified by substituting the total Ap/p for y(Thot - 
Tmeit). The density profile along the axis is given by 


Ap(z) 

P 


T(T m -T mai )e- Z/6 ’- + .16 (P"^ ggsl fTzT 


^-z/8 


PHg M Zn 18 Te 


where it is assumed that the density of the alloy is the weighted average of the 
two components. The composition profile falls to its half-value at 1.03 cm, thus X 
is taken to be 0.56 cm. 


From these values, u max = 31.63 (g/go) cm/sec, C = 1.28x10 4 (g/g 0 )C„ / k, and £ 

= 2.6x1 0 4 (g/go)- Therefore, for £ < 1%, the transverse acceleration must be 
<0.38 micro-g. 


4.4. Protein Crystal Growth 

Protein crystals grown from solution will be accompanied by a convective growth 
plume that flows along the side of the growing crystal as solute is taken up 
leaving a lower density depletion boundary layer. The velocity of the plume and 
the width of the boundary layer (the distance away from the surface where the 
velocity is maximum) can be estimated [2] from the approximate solution of 
Goldstien to the classical problem of free convective flow along a vertical heated 
plate [3]. 
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v max = 0.776 
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8 = 1.31 
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V 
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where h is the height of the crystal, Sc is the Schmidt number (v/D), and the 
Grashoff number Gr =g(Ap/p)h 3 / v 2 . For a large protein crystal assume 

h =0.1 cm 
v = 0.01 cm/sec 2 
D = 1 0-6 cm/sec 2 
Ap/p = 0.01. 


The Gr = 98 (g/go), the boundary layer thickness is 8 = 42(g/go)' 1/4 pm and v max = 
77 (g/go) 1 72 pm/sec. In order to have diffusion controlled transport, this boundary 

layer thickness must be » the chemical diffusion length h/2 = 500 pm. This 
requires (g/g 0 ) « (42/500) 4 = 5x1 O’ 5 . Therefore, in order to assure diffusive 
controlled growth of protein crystals, the continuous residual acceleration 
must be < 5 micro-g’s. 

The constant attitude required for the CGF experiments may have an unwanted 
effect on protein crystal growth in that the crystals will continue to migrate in the 
same direction at their Stokes velocity given by 

_ 2 g(Ap/p)a 2 

v Stokes gv 

For the above case, assuming a density difference between the crystal and the 
solution of 0.1 gm/cm 3 , vs t0 kes = 5 (g/go) cm/sec. Even at a steady acceleration 
of 1 micro-g, the crystals will migrate 1 cm in slightly more than two days. 
This will cause them to grow against the container walls which is believed 
to be undesirable from the stand point of growth conditions. This has not 
been a problem on previous missions because the attitude was not held 
constant for long enough periods of time for this to occur. 

Since the response time to a diffusion process goes as a 2 /v for momentum 
diffusion and a 2 /D for chemical diffusion, protein crystal growth may be more 
sensitive to the higher frequency accelerations because of the much shorter 
lengths scales involved. For growing crystals that are 100 pm in diameter, the 
response time for flow is only 2.5 milliseconds. Therefore, periodic disturbances 
with frequencies less than 64 Hz will produce flows approaching their steady 
state magnitude in the vicinity of the growing crystal. The chemical diffusion 
times are of course much slower, and the effects may not be felt with regard to 
diffusive transport. However, if flows of a few 10’s of microns/sec actually 
affect the attachment kinetics as some experiments have suggested [4,5], 


137 



milli-g accelerations with frequencies in the 1-60 Hz range could have an 
effect on the growth of such crystals. 


4.5 Marangoni Convection in Closed Containers 

The Marangoni Convection in Closed Containers experiment consists of a 
partially filled cylindrical tube (length L = diameter) containing water or a 
fluorocarbon oil. A longitudinal thermal gradient is applied an the behavior of the 
fluid is to be observed. The buoyancy contribution to any flows that might be 
observed is given by 

= g(Ap/p)L 2 
max 144V3 

For L = 2.54 cm and AT = 20°C, the Ap/p for water is 0.0068. The u ma x is 17.32 
(g/g 0 ). At 1 micro-g the flow will be 0.17 microns/sec. In an observation time of 5 
minutes, the marker particles will move only 52 microns - approximately half their 
diameter. Therefore, the contribution from buoyancy driven flows will be 
undetectable in this experiment. 0 
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SECTION 5 


RESULTS AND RECOMMENDATIONS 

5.1 USML-1 RESULTS 

As stated previously, the estimates of the g-levels required to effectively eliminate 
radial segregation in the furnace experiments that were made prior to the launch 
of USML-1 were based of the best model available at that time. Since it was 
expected that the USML-1 would be flown such that the residual acceleration 
would be along the furnace axis, the major attention up to that point had been 
given to predicting radial segregation caused by accelerations in the axial 
direction. In the post mission analysis, it became clear that unexpected forces, 
thought to be associated with the operation of the flash evaporator system, 
produced significant transverse accelerations. This was first noticed when the 
direction of the ball motion in the passive accelerometer (PAS) did not coincide 
with the furnace axis [3]. The OARE accelerometer indicated a transverse 
acceleration of 0.5 - 1 .0 micro-g which turned out to be somewhat larger than the 
axial acceleration due primarily to atmospheric drag [2]. This prompted the 
development of a more accurate model for the effects of transverse accelerations 
on Bridgman growth experiments which is described in Section 3. Using the 
results from this model, a new set of predictions for the melt growth experiments 
for calculated and are displayed in Table 5.1 . As can be seen by comparing the 
prediction in Section 4 with Table 5.1, the more accurate model requires more 
stringent control of the transverse accelerations, particularly for the case of large 
Sc as in the case of Lehoczky’s experiment. 

5.1.1 Lehoczky’s Experiment 

Lehoczky’s experiment was prematurely terminated shortly after the time at 
which steady state growth was achieved. However, the portion of the sample 
grown under steady state conditions strongly suggested the presence of an 
unanticipated transverse acceleration based on the interface shape, radial 
composition, and quenched-in dendrite structure. [3]. The radial segregation 
was measured to be 0.4 on one side and 0.18 on the other. The predicted value 
was 0.192 is of the right order of magnitude, but is somewhat less than 
observed. The primary difference is believed to be due to the fact that the model 
assumes a flat interface perpendicular to the growth ampoule; whereas, in the 
actual experiment, the interface was tilted due to the tendency of the denser 
rejected solute to slump to the side of the ampoule in the direction of the 
transverse acceleration. 

5.1.2 Larson’s Experiment 

Larson’s experiment is relatively insensitive to convective mixing because the k- 
value is close to 1 and, in fact, it is possible to achieve near diffusion control of 
solute distribution on the ground. Also, the melt did not completely fill the 
ampoule and, because of the unexpected transverse acceleration, the melt 
solidified along one side of the ampoule leaving a free surface along the other 
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Table 5.1 


Sensitivity of USML-1 experiment to 1 micro-g transverse acceleration. 



Matthiesen 

Matthiesen 

Larson 

Lehoczky 


Se:GaAs 

Se:GaAs 

Cd.95Zn.05Te 

Hg. 84 Zn.i 6 Te 


0.75 

0.75 

0.75 

0.4 

Ap/p 

0.0085 

0.0085 

0.01 

0.0956 

A. (cm) 

1.7 

1.7 

2.0 

1.6 

k 

0.1 

0.1 

1.2 

4.35 

v (cm 2 /sec) 

0.003 

0.003 

0.00435 

0.008 

RlBllillBi: 

10- 4 

10-4 

10-4 

6 x 10- 6 

vg (cm/sec) 

2 . 5 x 10- 4 

5x10-4 

4.4 xIO " 5 

4.1 xIO ' 6 

Gr 

0.39 

0.39 

0.218 

0.094 

Sc 

30 

30 

43.5 

1333 

8 (cm) 

0.4 

0.2 

2.273 

1.463 

a/X 

0.441 

0.441 

0.375 

0.249 

a/5 

1.875 

3.75 

0.33 

0.273 

umax (cm/sec) 

1.9x10 - 5 

1.9x10-5 

1.45 xIO ' 5 

1.75 xIO ' 5 


1.5x10 - 5 

1.5x10-5 

1.1 xIO 5 

1.23 xIO 5 

c 

0.047 

0.042 

0.002 

0.192 
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side. This violates the thermal and fluid no-slip boundary conditions assumed in 
the model and allows the presence of Marangoni flow which could seriously 
affect the solute redistribution. Mention was made of the fact that radial 
segregation was disturbed by the asymmetric acceleration and thermal fields in 
the flight sample, but details have not yet been reported [4], 

5.1.3 Matthiesen’s Experiment 

Matthiesen used a spring loaded PBN plunger to eliminate free surfaces in his 
experiment [5], The axial dopant distribution in his first ingot indicated diffusion- 
controlled growth initially, but indicated complete mixing after the first rate 
change. The second ingot exhibited complete mixing throughout the entire 
growth regime. In both cases, large voids were found in the solid, apparently 
produced by gas bubbles in the melt. Radial segregation was reported to be 
from 0.3-0.5 in both the ground control as well as in the flight samples. This 
appeared to be primarily due to the large interface curvature imposed by the 
thermal profile. 

It is difficult to understand how transverse accelerations of the order of 1 micro-g 
could produce flows greater than the translation velocity which would be required 
for complete mixing to occur. It is more likely that temperature gradients across 
the bubbles that formed produced Marangoni flows that could easily account for 
the complete mixing of the solute [6). 

5.1.4 Other experiments 

None of the other primary experiments on USML-1 offer the possibility of testing 
the models developed in this study. Vapor crystal growth experiments, such as 
performed by Wiedemeier, are relative immune to low level accelerations 
because of the high kinematic viscosity and low Schmidt numbers of vapors 
which greatly reduces the Gr Sc product. Estimates were made of the 
accelerations required to achieve diffusion controlled growth conditions in the 
protein crystal growth experiments of DeLucas et al. These were apparently 
achieved and some excellent crystals were grown, but since the influence of low 
level accelerations on the perfection of protein crystals growing from solution is 
not understood, these experiments offer no definitive test of the models. Similar 
statements apply to the zeolite crystal growth experiments. 

The surface tension driven convection experiment of Ostrach was so completely 
dominated by surface tension forces that it would be virtually impossible to 
extract any buoyancy-driven effects. Small surface fluctuations were reported 
that appeared to correlate with the vernier thruster firings, but these did not affect 
the observed flows [7]. 

The drop physics experiments only required g-levels low enough to be 
neutralized by the acoustic positioning field which apparently was the case on 
USML-1. 
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For the most part, the Glovebox experiments were qualitative rather than 
quantitative in nature. The Passive Accelerometer System (PAS) did provide 
independent verification of the direction and magnitude of the quasi-steady 
residual acceleration using a model based on Stokes’ law modified to account for 
wall effects [1], 

An attempt to observe flows resulting from transient accelerations from thruster 
firings was made using the cell containing water in the Marangoni Convection in 
Closed Containers (MCCC) experiment [6]. (Water surfaces are so easily 
contaminated with surface active impurities that the Marangoni effect is not 
observed unless extreme precautions are taken.) This cell, which contained 
marker particles illuminated by a light sheet was to have been placed in the back 
of the glove box and observed with a video camera for a number of hours while 
other experiments were being performed in the glove box. Unfortunately, this 
data was apparently not recorded for reasons that are not clear. 

Small crystals that were nucleated during the Nucleation of Crystals from 
Solution (NCS) experiments were observed to drift toward the chamber wall 
under the influence of microgravity [8], but no quantitative data were recorded. 

Several abrupt shifts in the particle distribution were seen during the long (20 
minute) runs of the Particle Dispersion Experiment (PDE) which were attributed 
to vernier thruster firings. These shifts did not affect the experiment, and, in fact, 
provided information in distinguishing which particles were free floating from 
those that were stuck to the wall [9]. 


5.2 SUMMARY 

It has been shown for dilute systems grown by the Bridgman method that the 
radial segregation produced either by an axial and by a transverse low-level 

quasi-steady acceleration scales as (1-k) Gr Sc f(a/A,, a/8, and k) where Gr is 
calculated with the radial thermal gradient in the case of axial acceleration, and 
with the axial thermal gradient in the case of transverse acceleration. (The 

function of (aA., a/8, and k) is of course different for the two cases.) Since the 
axial gradients are generally greater than the radial gradient, transverse 
accelerations are considered more critical. The transverse result also applies to 
the case of a non-dilute system provided the Gr Sc is small. 

For the case of a dilute system with a linear axial thermal gradient, Garandet [10] 
showed that the radial segregation caused by a steady transverse acceleration is 
given by 


c = 


(1-k)GrScPe g 

1766 


Pe g « 1 


(5.1) 


and 
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, P©g >>1 


(5.2) 




(l-k)GrSc 
k 4 


Pe„ 


where Gr = g p (AT/L) W 4 / v 2 and Pe g = v g W / D =2 a / 8. 

Similar results were obtained in Section 2 and it was also show that the effect of 
a gradient that falls off exponentially can produce several times the radial 
segregation as a linear gradient with the same temperature difference. For 
maximum coupling between the flow field and the diffusion field, it was found that 


c 


GrSc 

1428 


(1-k), 


k < 1 


and 


= GrSc(1-k) k>1 
1654 k 


(5.3) 


(5.4) 


where here Gr = g (3 AT W 3 / v 2 . 

Comparing this to Eq. (1 .5), we see that the Langbein and Tiby model over- 
predicts the variation in solute redistribution by more than three orders of 
magnitude if their length scale is taken as the ampoule width. 

Even so, it may be seen that achieving diffusion limited transport in growth from 
the melt will be difficult since the Gr Sc must be held to 0(10). For a dilute 

system such as Se:GaAs, 8p/p = 0.0085, v = 0.003 cm 2 /sec, Sc = 30, and k = 

0. 1 . If W = 1 .5 cm, Gr = 31 87 g and £ = 67 g. If the radial segregation is to be 
held to < 1%, g < 0.15 micro-g. This can be reduced somewhat by detuning the 
experiment away from the conditions of maximum coupling (as may be seen by 
comparing this result to the estimates for Matthiesen’s experiment in Table 4.1), 
but the range over which this is possible is limited by practical considerations. 

On one hand, one wants the v g as large as possible to minimize the growth 
transient and to be able to grow as much as possible in the time available. 
However, increasing the growth rate increases radial segregation as long as P g < 

1 . One could lower the sensitivity to radial segregation by increasing Pe g »1 to 

get on the back side of the curve where £ ~ Pe g - 3 ; however, this is constrained by 
thermal gradient if one wants to prevent interfacial breakdown from constitutional 
undercooling. Using Tiller’s constitutional supercooling criteria, this constrains 
the growth Peclet number by 

_ v gW _W WVT (1-k) _ W(AT /Jk) (1-k) (5 5) 

P0g D 8 < mC„ k mC„ k 

where m is the slope of the liquidus. 
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Attempting to increase the allowable growth rate by increasing the thermal 
gradient increases the driving force for segregation, causes difficulty in controlling 
the interface shape, and produces more thermal stress in the newly solidified 
materials which promotes defect formation. Better quality crystals can be 
obtained by increasing X and lowering the growth rate, but this increases the time 
and ampoule length required to growth a reasonable amount of material under 
steady state growth conditions. 

The situation is even more demanding for the non-dilute ll-VI alloy systems 
which can have Sc = 0(1 0 3 ). Here the X = 8 since the density gradient is 
determined by the diffusion field. Extremely slow growth rates are required 
because of the small values of D that are characteristic of some of these 
systems. Again, to grow a reasonable amount of material in the time available, 

the thermal gradient must be adjusted so that Pe g = W/8 = 0(1 ). For Pe g = 1 , the 
segregation can be expressed as 


GrSc(l-k) 
^ 3200 


k < 1 


and 


GrSc (k-1) 
1778 k 


k> 1. 


(5.6) 


(5.7) 


For Ap/p = 0.0956, Sc = 1333, and v = 0.008 cm 2 /sec, k = 4.35, and W = 0.8 cm; 
Gr = 764 g and £ = 441 g. Now to attain < 1% segregation, the residual 
acceleration along the transverse axis would have to be 23 nano-g’s! Again 
some relief can be obtained by reducing the growth rate to lower the Pe g as may 
be seen in Fig. 3.12. The Pe g in Lehoczky’s experiment was 0.547 which 
accounts for the lower value in Table 5.1 . However, to be able to grow any 
appreciable amount of material in a reasonable length of time, even using a pre- 
profiled solute distribution, this growth rate is approaching the lower practical 
limit. 


5.3 RECOMMENDATIONS 

For future low gravity experiments, it would be desirable to obtain some relief to 
this very stringent requirement for limiting the transverse accelerations. If this 
requirement applied to crystals grown on earth by the vertical Bridgman process, 
it would be necessary to maintain the furnace axis to better than 27 nanoradians 
in order to avoid assymmetric radial segregation. Obviously, the vertical density 
gradient seeks the local vertical which eliminates the tangential component. In 
the absence of any axial acceleration, which was assumed here, there is no 
stabilizing axial gradient. 
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5.3.1 Density Gradient Stabilization 

Therefore, one possibility for reducing this transverse acceleration requirement is 
to purposely add a stabilizing axial acceleration. This possibility was analyzed in 
Section 1.6 where it was shown that the effect of an axial acceleration was to 
alter the maximum convective by the factor IRa c l/(IRa c l+Ra(x)). If we define 
Ra(x) as g x (Ap/p) W 4 / (v k L), then IRa c l = 500.56. Now for dilute systems, the 
segregation is proportional to 


£ - GrSc 


500 ] 

f 

~ GrSc 

V 500 + Ra(x) J 

V 


500 

500 + (g x /g y )(W/\)GrPr 


(5.8) 


Since Gr is small in a microgravity experiment, it is obvious that there is little 
benefit to be gained by applying an axial acceleration to stabilize the flow for 
systems with low Pr. However, for non-dilute systems in which the density 
gradient is predominantly due to solutal effects, the Pr in the above equation is 
replaced by Sc. For a non-dilute systems with k > 1 the maximum radial 
segregation from transverse accelerations can be written from Eq. (5.6) as 


GrSc(k-l), 

r 500 ^ 

| GrSc(k-l) 

' 500 ^ 

1778 k ' 

V 500+Ra(x) j 

' 1778 k 

[ 500 + (g x / g y )( W / 8) Gr Sc J 


The stabilizing effects of axial acceleration at different gravity levels for 
Lehoczky’s experiment can be seen in Figure 5.1 . If the furnace axis is oriented 
within ±5° of the residual g-vector for accelerations < 10‘ 7 go, the primary effect 
will be just the reduction of the transverse acceleration by sin 0. Increasing the 
acceleration above this value increases the segregation, but at a much slower 
rate due to the stabilizing effect of the axial gradient. It is interesting to note that 
by changing the geometry in such a way as to increase both the Gr and Ra, (e.g. 
by doubling the width of the ampoule, one can actually reduce the radial 
segregation caused by buoyancy-driven flows resulting from axial density 
gradients at the higher g-levels. However, as the axial acceleration is increased, 
flows driven by radial density gradients will also increase and possibly negate this 
effect. 
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Fig. 5.1 Radial segregation (zeta) as a function of the magnitude of the residual 
acceleration vector for different orientations relative to the furnace axis (theta) 
and different ampoule widths (W). For W = 0.8 cm used in Lehoczky’s 
experiment, orienting the furnace axis to within ± 5° of the residual acceleration 
produces approximately an order of magnitude reduction in segregation at 0.1 
micro-g primarily because of the sin (theta) reduction of the transverse 
component of the acceleration. The reduction becomes more pronounced as the 
total acceleration is increased as the stabilizing effect of the axial gradient 
becomes more important; however, the segregation continues to increase with 
increasing rate, although at a lessor rate. Doubling the width of the ampoule 
increases both the Gr and the Ra which makes the stabilizing effect of the 
density gradient more effective at a lower acceleration and can actually result in 
less segregation as can be seen by the cross-over of the two curves at ~ 1 micro- 

g- 
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5.3.2 Magnetic Stabilization 

The other possible way to stabilize such a system is the application of magnetic 
fields. This is illustrated in Fig. 5.2 which shows the maximum flow velocity in a 
vertical cell with width W and height 2 W with differentially heated side walls in 

unit gravity and in microgravity. The AT was taken as 10° C to roughly represent 
the radial gradient in a typical semiconductor melt using the vertical Bridgman 
system. The dotted diagonal lines represent Re = 2000, above which unsteady 
or turbulent flow develops; Peihermai = 1 . which marks the division between 
conductive and convective dominated heat transfer, and Pe SO iutai = 1, which 
marks the division between conductive and convective dominated mass transfer. 

The horizontal branches are the results of the application of a vertical magnetic 
field of strength B (Tesla) to a melt whose conductivity is a (Siemens or mho/m). 
The value a B 2 = 3.6 x 10 6 represents an approximate upper limit to magnetic 

stabilization (o = 10 5 S, B = 6 T) of semiconductor melts. With this amount of 
damping it is possible to avoid unsteady flows in even the largest melts in unit 
gravity and maintain conductive heat transfer, but convective mass transfer will 
still predominate unless the width is reduced to micron dimensions. 

In microgravity, the magnetic damping is much more effective due to the lower 
driving force. However, it may be seen that for melts with o = 0(1 0 5 S), a field of 
> 0.1 T (1 KG) would be required to obtain a significant reduction in flows for 
widths = 0(1 cm). 

It was shown in Chapter 1 .7 that if a strong magnetic field applied along the axis 
of a horizontal slot, the maximum flow along the walls was given by 


. _ Gr(W/L)( L V s v 

Umaa " 6^3 Ha 2 lwJ W’ 


Ha »1. 


(5.10) 


The ratio of this to the maximum velocity in a horizontal slot with a gradient AT/L 
is 


^ = )■ Ha » 1. (5.11) 

u Ha IwJ 

From this, it may be seen that to obtain a given reduction in flow from magnetic 
damping, the Ha - L/W. The reason for this behavior may be understood by 
examining Fig. 1 .23. Since the field inhibits flow in the transverse direction, the 
flow field adjusts itself so that the transverse flow is spread as much as possible 
throughout the ampoule to minimize the transverse velocity. Thus lengthening 
the ampoule gives the flow more room to adjust and allows a higher maximum 
axial flow velocity. This does have an added beneficial effect in that the 
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Chamber width (cm) 


Fig. 5.2 Map of maximum velocity in a vertical chamber with a 2:1 aspect ratio 
and differentially heated side walls. The Pr is taken to be 0.01 , v = 0.0013 
cm 2 /sec, (3 AT= 2.5x1 0 3 . Also shown are lines representing Re = 2000, above 
which flow is expected to be turbulent; Pe(Ther) = 1 , which marks the division 
between conductive and convective heat transport; and Pe(sol) = 1 , which marks 
the division between diffusive and convective mass transport for Sc = 10. 
Horizontal lines indicate the velocity of a fluid with conductivity cr (Siemens) in the 
presence of a vertical magnetic field B (Tesla). 
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maximum axial flow is moved further from the growth interface, thus reducing the 
coupling between the flow field and the concentration gradient which will result in 
less solute redistribution. 

Scaling from Table 5.1 , we see that in order to attain £ = 0.01 in Lehoczky’s 
experiment with 1 micro-g transverse acceleration, a reduction of 20 in Gr will be 
required. Using Eq. 5.1 1 to obtain a rough estimate of the required Ha, we find 
that Ha = 15.5 LAN. For LAN =5, p = 5000 Kg/m 3 , v = 8x1 O' 7 m 2 /sec, a = 0.004 
m, the a B 2 product would have to be 6x1 0 4 {LAN) 2 m.k.s. units. Estimating the 
conductivity of the melt to be 10 5 S, the required B is 0.775 (L/W)T or 38.75 KG. 
This may be an over-estimate because it assumes a linear gradient over the 
length of the melt and does not account for the reduced coupling between the 
flow field and the composition gradient near the interface. Therefore, a more 
detailed calculation is required to give a more definitive estimate. However, it 
seems clear that a very large magnetic field will be needed, probably requiring a 
superconducting magnet. 

If the orientation of the transverse acceleration is known a priori, it may be 
possible to reduce the requirement on B by the ratio (L/W) by using as transverse 
magnetic field with B oriented perpendicularly to the acceleration vector since it is 
only necessary to control the flow in one plane. However, as may be seen in Fig. 
1 .26, this has the result holding the axial flow close to the walls of the ampoule 
until the region near the growth interface. Thus, some of the advantage may be 
lost by the increased coupling between the flow field and the composition 
gradient. 

Magnetic damping coupled with microgravity would appear to have great promise 
in being able to extend diffusion control of less demanding (lower Sc) systems to 
larger diameters since Ha ~ width and the convective velocity becomes 
insensitive to size scale for Ha » 1 . The segregation behavior of a system 
whose density gradient falls exponentially in an axial magnetic field needs to 
analyzed in more detail such as was done in Section 3 in order to asses the 
effects of the reduced coupling between the flow field and the composition 
gradient. 

5.4 CONCLUSIONS 

Approximate analytical models for the flow fields in various geometries have 
been developed and verified against computational results. These allow the 
accurate prediction of velocities in a wide variety of microgravity experiments as 
well as in ground based experiments as long as the Reynolds numbers are low 
enough so that inertial terms can be neglected. Perturbation models for solute 
redistribution in Bridgman melt growth systems have been coupled to these flow 
models which allow the prediction of radial segregation due to axial as well as 
transverse accelerations. Again these have been verified by testing against 
computational results. 
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Unfortunately, experimental difficulties with the various crystal growth 
experiments using the Bridgman-type Crystal Growth Furnace (CF) did not allow 
a definitive experimental test of these models, but it was obvious that these 
experiments were affected by the unexpected transverse acceleration that 
ranged between 0.5 and 1 micro-g. 

The ability to achieve diffusion controlled solute redistribution in non-dilute 
systems with large Sc remains a formidable challenge, even in a microgravity 
environment. Efforts to align the furnace axis with the local acceleration vector 
certainly help in reducing the transverse component of the acceleration, but 
effects such as diurnal variations in the atmosphere as well as uncontrolled leaks 
and vents place practical limitations on how well this may be accomplished. 
Density gradient stabilization the application of magnetic fields may help 
overcome this problem, but more research will be required to determine their 
effectiveness. 
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